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Chapter 1 
Formation channels of single and binary 
stellar-mass black holes 


Michela Mapelli 


Abstract These are exciting times for binary black hole (BBH) research. LIGO and 
Virgo detections are progressively drawing a spectacular fresco of BBH masses, 
spins and merger rates. In this review, we discuss the main formation channels of 
BBHs from stellar evolution and dynamics. Uncertainties on massive star evolution 
(e.g., stellar winds, rotation, overshooting and nuclear reaction rates), core-collapse 
supernovae and pair instability still hamper our comprehension of the mass spectrum 
and spin distribution of black holes (BHs), but substantial progress has been done 
in the field over the last few years. On top of this, the efficiency of mass transfer 
in a binary system and the physics of common envelope substantially affect the 
final BBH demography. Dynamical processes in dense stellar systems can trigger 
the formation of BHs in the mass gap and intermediate-mass BHs via hierarchical 
BH mergers and via multiple stellar collisions. Finally, we discuss the importance 
of reconstructing the cosmic evolution of BBHs. 


Key words: black holes — gravitational waves — pair instability — massive stars 
— binary systems — binary black holes — intermediate-mass black holes — stellar 
dynamics — star clusters — merger rates 


1.1 Introduction: Observational facts about gravitational waves 


On 2015 September 14, the LIGO interferometers captured the gravitational 
wave (GW) signal from a binary black hole (BBH) merger [5]. This event, named 
GW150914, is the first direct detection of GWs, about hundred years after Ein- 
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stein’s prediction. Over the last five years, LIGO and Virgo witnessed a rapidly 
growing number of GW events: the second gravitational wave transient catalogue 
(GWTC-2, [16]) consists of 50 binary compact object mergers from the 
first (O1), the second (O2) and the first part of the third observing run (O3a) of the 
LIGO-Virgo collaboration (LVC). Based on the results of independent pipelines, 
[383], [350], and claimed several additional GW candidates from O1 
and O2. Furthermore, several dozens of public triggers from the second part of the 
third observing run (O3b) can be retrieved at/https: //gracedb.ligo.org/| 

Among the aforementioned detections, GW170817, the first binary neutron star 
(BNS) merger detected during O2, is the first and the only GW event unquestionably 
associated with an electromagnetic counterpart to date 
: 

Several other exceptional events were observed during O3a: the first unequal- 
mass BBH, GW 190412 [17]; ii) the second BNS, GW190425 [6]; iii) the first black 
hole (BH) — neutron star (NS) candidate, GW190814 [19], and iv) GW190521, 
which is the most massive system ever observed with GWs [181[20]. Also, GW190521 
is the first BBH event with a possible electromagnetic counterpart [157]. This grow- 
ing sample represents a "Rosetta stone" to unravel the formation of binary compact 
objects. 


Astrophysicists have learned several revolutionary concepts about compact ob- 
jects from GW detections. Firstly, GW 150914 has confirmed the existence of BBHs, 
i.e. binary systems composed of two BHs. BBHs have been predicted a long time 
ago (e.g. [53]), but GW150914 is their first 
observational confirmation [2]. Secondly, GW detections show that a number of 
BBHs are able to merge within a Hubble time. Thirdly, GW170817 has confirmed 
the connection between short gamma-ray bursts, kilonovae and BNS mergers [11]. 

Finally, most of the LIGO-Virgo BHs observed so far host BHs with mass in 
excess of 20 Mo. The very first detection, GW 150914, has component masses equal 
to m, = 35.6'31 Mo and m = 30.62? Mo [3]. This was a genuine surprise for 
the astrophysicists [2], because the only stellar BHs for which we have a dynamical 
mass measurement, i.e. about a dozen of BHs in X-ray binaries, have mass < 20 Mo 
(2641[266] 249]. Moreover, most theoretical models available five years ago did not 
predict the existence of BHs with mass mgu > 30 Mo (but see 
for a few exceptions). Thus, the first GW detections have 
urged the astrophysical community to deeply revise the models of BH formation 
and evolution. 

The recently published O3a events add complexity to this puzzle. In particular, 
the secondary component of GW190814, with a mass m = 2501015 Mo [19], 
is either the lightest BH or the most massive NS ever observed, questioning the 
proposed existence of a mass gap between 2 and 5 Mo ([266][116). 

The merger product of GW190521 (m; = 1427?* Mo, [20]) is the first 
intermediate-mass BH (IMBH) ever detected by LIGO and Virgo and is the first 
BH with mass in the range ~ 100 — 1000 Mo ever observed not only with GWs but 
also in the electromagnetic spectrum. Moreover, the mass of the primary compo- 
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nent, m, = Caen Mo, falls inside the predicted pair instability mass gap, as we will 
discuss later in this review. According to the re-analysis of (see also (120]), 
GW190521 could be an intermediate-mass ratio inspiral with primary (secondary) 
mass 168*5 Mo (162 Mo), when assuming a uniform in mass-ratio prior. This 
interpretation avoids a violation of the mass gap, but requires the formation of an 
IMBH to explain the primary component. 


LIGO and Virgo do not observe only masses, they also allow us to extract in- 
formation on spins and merger rates. The local BBH merger rate density inferred 
from GWTC-2 is 23.917? (58735) Gpe~> yr-! within the 90% credible inter- 
val, when GW190814 is included in (excluded from) the sample of BBHs [15]. 
The LVC also estimated an upper limit for the merger rate density of BH-NS bi- 
naries (Žguns < 610 Gpc ? vit. (8]) and inferred a BNS merger rate density 
ÜfgNs = 32000 Gpc ^? yr^! within the 90% credible interval [15]. 

The uncertainties on the spins of each binary component are still too large to 
draw strong conclusions, apart from very few cases (for example, the spin of the 
primary component of GW190814 is close to zero, see Figure 6 of [19]). However, 
LIGO and Virgo allow to give an estimate of two spin combinations, which are 
called effective spin (Yer) and precessing spin (Xp). The effective spin is defined as 


ya = mZ +m x2) L (1.1) 


mı +m L' 


where m and m (y | and %2) are the masses (dimensionless spin parameters) of the 
primary and secondary component of the binary, respectively and L is the Newto- 
nian orbital angular momentum vector of the binary. The dimensionless spin param- 
eters are defined as 7; = S;c/(G m2), where S; is the magnitude of the spin vector, 
c the speed of light and G the gravity constant. Hence, the effective spin can take 
any values between —1 and 1, where y.s = 1 (—1) means that the two BHs are both 
maximally rotating and perfectly aligned (anti-aligned) with respect to the orbital 
angular momentum of the BBH, while y.; = 0 indicates either that the two BHs are 
perfectly non-rotating or that both spins lie in the plane of the BBH orbit. 
The precessing spin is defined as 


Xp max (B1 S11, B2821), (1.2) 


7 Bı Gm; 
where B, = 2+3q/2 and B? = 2+3/(2q), q — m/m; < 1, $4, and $5, are the 
components of the spin vectors perpendicular to the orbital angular momentum. 
Hence, x, can take values between 0 (no spin components in the orbital plane) and 
1 (at least one spin being maximal and lying in the orbital plane). %p is called pre- 
cessing spin because spin components misaligned with respect to the orbital angular 
momentum of the binary drive precession [313]. 

In current GW observations, small values of Xeff are preferred and %p is uncon- 
strained [7], with a few exceptions. GW151226, GW170729, GW190412 and 
GW190425, together with a few candidate events (e.g., GW190517.055101) show 
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support for positive values of yer [8][17][14]. If we look at the overall population of 
GWTC-2 BBHs [15], ~ 12% to 44% of BBH systems have spins tilted by more than 
90? with respect to their orbital angular momentum, supporting a negative effective 
spin parameter. 

The precessing spin of GW190814 has a strong upper bound x < 0.07 [19]. 
Finally, GW190521 shows mild evidence for a non-zero precessing spin (yp = 
0.68+9:35 within the 90% credible interval, [20]). In contrast, spin measurements in 
X —ray binaries point to a range of spin magnitudes, including high spins (248][249]. 

This review discusses the formation channels of BHs and BBHs in light of the 
challenges posed by recent GW detections. This field has witnessed an exponential 
growth of publications and models in the last few years. While I will try to give an 
overview as complete as possible of the main models and connected issues, it would 
be impractical to mention every interesting study in the span of this review. 


1.2 The formation of compact remnants from single stellar 
evolution and supernova explosions 


Black holes (BHs) and neutron stars (NSs) are expected to form as remnants of 
massive ( Z 8 Mc) stars. An alternative theory predicts that BHs can also form 
from gravitational collapse in the early Universe (the so called primordial BHs, e.g. 
(721[58] [711 [173]. In this review, we will focus on BHs of stellar origin. 

The mass function of BHs is highly uncertain, because it may be affected by a 
number of barely understood processes. In particular, stellar winds and supernova 
(SN) explosions both play a major role on the formation of compact remnants. Pro- 
cesses occurring in close binary systems (e.g. mass transfer and common envelope) 
are a further complication and will be discussed in the next Section. 


1.2.1 Stellar winds and stellar evolution 


Stellar winds are outflows of gas from the atmosphere of a star. In cold stars (e.g. 
red giants and asymptotic giant branch stars) they are mainly induced by radiation 
pressure on dust, which forms in the cold outer layers (e.g. [354]). In massive hot 
stars (O and B main sequence stars, luminous blue variables and Wolf-Rayet stars), 
stellar winds are powered by the coupling between the momentum of photons and 
that of metal ions present in the stellar photosphere. A large number of strong and 
weak resonant metal lines are responsible for this coupling (see e.g. fora 
review). 

Understanding stellar winds is tremendously important for the study of compact 
objects, because mass loss determines the pre-SN mass of a star (both its total mass 
and its core mass), which in turn affects the outcome of an SN explosion 
[225]|232] 48]. 
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Early work on stellar winds (e.g. [200]) highlighted that the mass loss 
of O and B stars depends on metallicity as ri œ Z“ (with œ ~ 0.5 — 1.0, depending 
on the model). However, such early work did not account for multiple scattering, 
i.e. for the possibility that a photon interacts several times before being absorbed 
or leaving the photosphere. Vink et al. (2001, (360]) accounted for multiple scat- 
terings and found a universal metallicity dependence ri? œ Z0? v^, where væ is the 
terminal velocity and p — —1.23 (p — —1.60) for stars with effective temperature 
Tett Z 25000 K (12000 K Z Tere £ 25000 K). 

The situation is more uncertain for post-main sequence stars. For Wolf-Rayet 
(WR) stars, i.e. naked helium cores, predict a similar trend with metallic- 
ity ri œ Z096. With a different numerical approach (which accounts also for wind 
clumping), find a strong dependence of WR mass loss on metallicity but also 
on the electron-scattering Eddington factor I; = K,L/(47t c Gm), where Ke is the 
cross section for electron scattering, L is the stellar luminosity, c is the speed of 
light, G is the gravity constant, and m is the stellar mass. The importance of I; has 
become increasingly clear in the last few years [155][361] [359], but only few stellar 
evolution models include this effect. 

For example, adopt a mass loss prescriptions for massive hot stars (O 
and B stars, luminous blue variables and WR stars) that scales as ri œ Z*, where 
a = 0.85 if I; < 2/3, a = 2.45 - 241; if 2/3 € I; € 1 and a = 0.05 if I; > 
1. This simple formula accounts for the fact that metallicity dependence tends to 
vanish when the star is close to be radiation pressure dominated, as clearly shown 
by figure 10 of [154]. Figure [1.1] shows the mass evolution of a star with zero-age 
main sequence (ZAMS) mass mzams = 90 Mo for seven different metallicities, as 
obtained with the SEVN code [322]. At the end of its life, a solar-metallicity star 
(here we assume Zo = 0.02) has lost more than 2/3 of its initial mass, while the 
most metal-poor star in the Figure (Z = 0.005 Zo) has retained almost all its initial 
mass. 

Other aspects of massive star evolution also affect the pre-SN mass of a star. For 
example, surface magnetic fields appear to strongly quench stellar winds by mag- 
netic confinement [141][275] [180]. In particular, show that a non-magnetic star 
model with metallicity ~ 0.1 Zo and a magnetic star model with solar metallicity 
and Alfvén radius RA ~ 4Rc undergo approximately the same mass loss according 
to this model. This effect cannot be neglected because surface magnetic fields are 
detected in ~ 10 per cent of the hot stars [366], but is currently not included in 
models of compact-object formation. 

Finally, rotation affects the evolution of a massive star in several ways (e.g. 
(2171761[202]203][77]). As a general rule of thumb, rotation increases the stellar lu- 
minosity. This implies that mass loss is generally enhanced if rotation is accounted 
for. On the other hand, rotation also induces chemical mixing, which leads to the 
formation of larger helium and carbon-oxygen cores. While enhanced mass loss im- 
plies smaller pre-SN masses, the formation of bigger cores has strong implications 
for the final fate of a massive star, as we discuss in the following Section. 
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Fig. 1.1 Evolution of stellar mass as a function of time for a star with ZAMS mass mzams = 90 Mo 
and seven different metallicities, ranging from 0.005 Zo up to Zo (we assumed Zo = 0.02). These 
curves were obtained with the SEVN population-synthesis code [322], adopting PARSEC stellar 
evolution tracks [75]. 


1.2.2 Core-collapse supernova (SN) or direct collapse 


Whether a star undergoes a successful core-collapse SN or a failed SN is the first 
key question to address, in order to assess the properties of the final compact object. 
A star undergoing a successful core-collapse SN explosion will leave a NS or a light 
BH, while stars that end their life with a failed SN will become rather massive BHs 
(> 20 Mo), because most (if not all) of the final mass of the star collapses into a BH 
directly. Addressing this question is a challenge for several reasons. 

The mechanisms triggering iron core-collapse SNe are still highly uncertain. 
The basic framework and open issues are the following. As the mass of the central 
degenerate core reaches the Chandrasekhar mass [73], the degeneracy pressure of 
relativistic electrons becomes insufficient to support it against collapse. Moreover, 
electrons are increasingly removed, because protons capture them producing neu- 
trons and neutrinos. This takes the core into a new state, where matter is essentially 
composed of neutrons, which support the core against collapse by their degeneracy 
pressure. To reach this new equilibrium, the core collapses from a radius of few 
thousand km down to a radius of few ten km in less than a second. The gravitational 
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energy gained from the collapse is W ~ 5 x 10° erg (mpys/1.4Mc;)? (10km/Rpns), 
where mpys and Rpys are the mass and radius of the proto-neutron star (PNS). 

The main problem is to explain how this gravitational energy can be -at least 
partially— transferred to the stellar envelope triggering the SN explosion [57]. 
Several mechanisms have been proposed, including rotationally-driven SNe and/or 
magnetically-driven SNe (see e.g. and references therein). The most com- 
monly investigated mechanism is the convective SN engine (see e.g. [133]). Accord- 
ing to this model, the collapsing core drives a bounce shock. For the SN explosion to 
occur, this shock must reverse the supersonic infall of matter from the outer layers 
of the star. Most of the energy in the shock consists in a flux of neutrinos. As soon as 
neutrinos are free to leak out (because the shock has become diffuse enough), their 
energy is lost and the shock stalls. The SN occurs only if the shock is revived by 
some mechanism. In the convective SN scenario, the region between the PNS sur- 
face and the shock stalling radius can become convectively unstable (e.g. because 
of a Rayleigh-Taylor instability). Such convective instability can convert the energy 
leaking out of the PNS in the form of neutrinos to kinetic energy pushing the con- 
vective region outward. If the convective region overcomes the ram pressure of the 
infalling material, the shock is revived and an explosion is launched. If not, the SN 
fails. 

While this is the general idea of the convective engine, fully self-consistent sim- 
ulations of core collapse with a state-of-the-art treatment of neutrino transport do 
not lead to explosions in spherical symmetry except for the lighter SN progenitors 
(< 10 Mo, (121) [110]). Simulations which do not require the assumption of spher- 
ical symmetry (i.e. run at least in 2D) appear to produce successful explosions from 
first principles for a larger range of progenitor masses (see e.g. [254]. How- 
ever, 2D and 3D simulations are still computationally challenging and cannot be 
used to make a study of the mass distribution of compact remnants. 

Thus, in order to study compact-object masses, SN explosions are artificially 
induced by injecting in the pre-SN model some amount of kinetic energy (kinetic 
bomb) or thermal energy (thermal bomb) at an arbitrary mass location. The evolu- 
tion of the shock is then followed by means of 1D hydrodynamical simulations with 
some relatively simplified treatment for neutrinos. This allows to simulate hundreds 
of stellar models. 

Following this approach, O'Connor & Ott (2011, (262]) propose a criterion to 
decide whether a SN is successful or not, based on the compactness parameter: 


m/Mo 
Sm = R(m)/1000km’ (3) 
where R(m) is the radius which encloses a given mass m. Usually, the compactness 
is defined for m = 2.5 Mọ (&25). measure the compactness at core bouncd!] 
in their simulations and find that the larger & 5 is, the shorter is the time to form a 
BH (as shown in their Figure 6). This means that stars with a larger value of &5 
are more likely to collapse to a BH without SN explosion. The work by Ugliano et 


1 show that £» 5 is not significantly different at core bounce or at the onset of collapse. 
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al. (2012, (351]) and Horiuchi et al. (2014, [171]) indicate that the best threshold 
between exploding and non-exploding models is &5 5 ~ 0.2. 

Ertl et al. (2016, ([110]) indicate that a single criterion (e.g. the compactness) 
cannot capture the complex physics of core-collapse SN explosions. They introduce 
a two-parameter criterion based on 


m(s — 4) 


M4 = ——— — 
4 Mo 


dm/Mo | (14) 
s=4 


ae ne ar at 
where Ma is the mass (at the onset of collapse) where the dimensionless entropy per 
baryon is s = 4, and u4 is the spatial derivative at the location of M4. This choice is 
motivated by the fact that, in their 1D simulations, the explosion sets shortly after 
M; has fallen through the shock and well before the shell enclosing M4 -- 0.3 Mo 
has collapsed. They show that exploding models can be distinguished from non- 
exploding models in the u4 versus M4 [4 plane (see their Figure 6) by a linear fit 


y(x) 2 ki x 4- Ko (1.5) 


where y(x) = Ha, x = M4 Us, and kı and kz are numerical coefficients which depend 
on the model (see Table 2 of [110]). The reason of this behaviour is that u4 scales 
with the rate of mass infall from the outer layers (thus the larger u4 is, the lower the 
chance of the SN to occur), while M4 u4 scales with the neutrino luminosity (thus 
the larger M4 L4 is, the higher the chance of a SN explosion). Finally, [110] find that 
fallback is quite inefficient (< 0.05 Mo) when the SN occur. 

The models proposed by O'Connor & Ott (2011, (262]) and Ertl et al. (2016, 
110], see also [I11]) are sometimes referred to as the “islands of 
explodability" scenario, because they predict a non-monotonic behaviour of SN ex- 
plosions with the stellar mass. This means, for example, that while a star with a mass 
m = 25 Mo and a star with a mass m = 29 Mo might end their life with a powerful 
SN explosion, another star with intermediate mass between these two (e.g. with a 
mass m = 27 Mo) is expected to directly collapse to a BH without SN explosion. 
Thus, these models predict the existence of islands of explodability, i.e. ranges of 
mass where a star is expected to explode, surrounded by mass intervals in which the 
star will end its life with a direct collapse. 

The models discussed so far depend on quantities (2.5, M4, H4) which can be 
evaluated no earlier than the onset of core collapse. Thus, stellar evolution models 
are required which integrate a massive star till the iron core has formed. This is 
prohibitive for most stellar evolution models (with few remarkable exceptions, e.g. 
FRANEC and MESA [272]. 

Fryer et al. (2012, [133]) propose a simplified approach (see also [132][134][51]). 
They suggest that the mass of the compact remnant depends mostly on two quan- 
tities: the carbon-oxygen core mass mco and the total final mass of the star mgn. 
In particular, mco determines whether the star will undergo a core-collapse SN or 
will collapse to a BH directly (namely, stars with mco > 11 Mo collapse to a BH 
directly), whereas mg; determines the amount of fallback on the PNS. In this formal- 
ism, the only free parameter is the time to launch the shock. The explosion energy 
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is significantly reduced if the shock is launched >> 250 ms after the onset of the 
collapse (delayed SN explosion) with respect to an explosion launched in the first 
~ 250 ms (rapid SN explosion, [133]). 

While this approach is quite simplified with respect to other prescriptions, 
and show that there is a strong correlation between the final carbon-oxygen 
mass and the compactness parameter $» 5 at the onset of collapse, regardless of the 
rotation velocity of the progenitor star (see figure 1 of (236]). Thus, we can conclude 
that the simplified models by can effectively describe the overall trend of a 
collapsing star, although they do not take into account several details of the stellar 
structure at the onset of collapse. 

Recently, propose an interesting alternative approach. They integrate a 
large grid of naked carbon-oxygen (CO) cores to the onset of core collapse and es- 
timate the explodability of each model with the compactness & 2.5 and with 
the parameter M4 [111]. Naked CO cores are faster and simpler to evolve than 
full stellar models (with hydrogen and helium) and are less sensitive to metallicity. 
They made available their grid of simulations for implementation into population- 
synthesis codes. Other works (e.g. [222]) highlight the stochasticity of 
the final direct collapse or core-collapse SN. Finally, even the most advanced for- 
malisms to derive the explodability of massive stars should be taken with a grain of 
salt, because of the complexity of the processes involved in core collapse SNe and 
because of the simplifications still included in the models [68]. 

Even if we were in the conditions to tell if a given star undergoes a failed SN 
instead of a successful SN, this would not mean we can automatically infer the final 
mass of the compact object. In the case of a failed SN, the main uncertainty on the 
final compact object mass is represented by the fate of the envelope [236]. In fact, 
the envelope of a massive giant star is rather loosely bound and even a small energy 
injection can unbind a fraction of it. 

Fernandez et al. (2018, (118]) show that a 0.1 — 0.5 Mo neutrino emission during 
the PNS phase causes a decrease in the gravitational mass of the core, resulting in 
an outward going pressure wave (sound pulse) that steepens into a shock as it travels 
out through the star. This might cause the ejection of a fraction of the loosely bound 
stellar envelope. According to [118], the ejected mass is a monotonically decreasing 
function of the envelope compactness (Figure 6 of [118]), defined as 


Mec/Mo 
Ree [Ro ' 


fam (1.6) 


where Mec and Rec are the total mass and radius of the star at the onset of core- 
collapse. With this formalism, the ejected mass is up to a few Mo for red super-giant 
stars, and 1 Mo for more compact stars like blue super-giant stars and WR stars. 

The possibility of a direct collapse is supported by observations. A survey 
conducted with the Large Binocular Telescope to find quietly disappearing stars 
reported evidence for the disappearance of a ~ 25 Mo red super-giant 
star [23]. In addition, surveys of SNe indicate a dearth of red super-giant progenitors 
with mass > 20 Mo associated with Type Ip SNe [182]. 
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1.2.3 Pair instability and the mass gap 


If the helium core of a star grows above ~ 30 Mo and the core temperature is 
2% 7 x 10° K at the end of carbon burning, the process of electron-positron pair pro- 
duction becomes effective. It removes photon pressure from the core producing a 
sudden contraction of the carbon-oxygen core, before the formation of an iron core 
[123] [42] 288] [375]. For mye > 135 Mo, the contraction cannot be reversed and the 
star collapses directly into a BH [375]. If 135 2 mye Z 64 Mo, the collapse trig- 
gers an explosive burning of heavier elements, especially oxygen and silicon. This 
leads to a pair instability SN (PISN): the star is completely disrupted, leaving no 
compact remnant [164]. For 64 Z mye Z 32 Mo, pair production induces a series of 
pulsations of the core (pulsational pair instability), which trigger an enhanced mass 
loss [375]. At the end of this instability phase, the star finds a new equilibrium and 
evolves towards core-collapse: a compact object with non-zero mass is produced, 
less massive than we would expect without pulsational mass loss. 

The main effect of (pulsational) pair instability is to open a gap in the mass 
spectrum of BHs between approximately ~ 507) Mo and ~ 120 Mo 
[321] [147] [145] 240] [239] [330) [355]. The large uncertainty on the edges of the mass 
gap is due to our poor understanding of the physics of massive stars. In particular, 
and integrate pure-He stars with MESA and show that 
the uncertainties on the PCa, y) 160 reaction rate are responsible for a change 
in the lower-edge of the mass gap of more than ~ 15 Mo. 

Mapelli et al. (2020, [236]) show that assuming that the hydrogen envelope col- 
lapses to the final BH can lead to an increase of the lower-edge of the mass gap from 
~ 40 Mo to ~ 65 Mo. Finally, combine the uncertainties on the "^C(a, y)! 6O 
reaction rate with the uncertainties on the collapse of the hydrogen envelope and find 
that the mass gap progressively reduces from ~ 80 — 150 Mọ (for a rate computed 
with the standard !*C(a, y)!6O reaction rate —1 o) to ~ 92 — 110 Mo (for a standard 
rate — 20) and even disappears (for a standard rate —3 0) as an effect of convection 
and envelope undershooting. These uncertainties leave open the possibility that the 
primary mass of GW190521 is the result of stellar evolution ([86][52] [340] [339]. 


1.2.4 The mass of compact remnants 


The previous Sections suggest that our knowledge of the compact remnant mass 
is hampered by severe uncertainties, connected with both stellar winds and SNe. 
Thus, models of the mass spectrum of compact remnants must be taken with a grain 
of salt. However, a few robust features can be drawn. 

Figure[I.2]is a simplified version of Figures 2 and 3 of Heger et al. (2003 [163]). 
The final mass of a star and the mass of the compact remnant are shown as a func- 
tion of the ZAMS mass. The left and the right-hand panels show the case of solar 
metallicity and metal-free stars, respectively. In the case of solar metallicity stars, 
the final mass of the star is much lower than the initial one, because stellar winds 
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Fig. 1.2 Final mass of a star (mg, blue lines) and mass of the compact remnant (Mrem, red lines) 
as a function of the ZAMS mass of the star. The thick black line marks the region where mg; = 
mzams. Left-hand panel: solar metallicity stars. Right-hand panel: metal-free stars. The red arrow 
on the left-hand panel is an upper limit for the remnant mass. Vertical thin black line in the right- 
hand panel: approximate separation between successful and failed SNe at Z = 0. This cartoon was 
inspired by Figures 2 and 3 of Heger et al. (2003 [163]). 


are extremely efficient. The mass of the compact remnant is also much lower than 
the final mass of the star because a core-collapse SN always takes place. 

In contrast, a metal-free star (i.e., a population III star) loses a negligible fraction 
of its mass by stellar winds (the blue and the black line in Figure[I.2] overlap). AS 
for the mass of the compact remnant, Figure[I.2|shows that there are two regimes: 
below a given threshold (~ 30 — 40 Mo) the SN explosion succeeds even at zero 
metallicity and the mass of the compact remnant is relatively small. Above this 
threshold, the mass of the star (in terms of both core mass and envelope mass) is 
sufficiently large that the SN fails. Most of the final stellar mass collapses to a BH, 
whose mass is significantly larger than in the case of a SN explosion. The only 
exception is represented by the pair instability window: single metal-free stars with 
ZAMS mass mzams ~ 140 — 260 Mo undergo a PISN and are completely destroyed, 
while single metal-free stars with mzams ~ 110 — 140 M; undergo pulsational pair 
instability and leave smaller compact objects. In this simplified cartoon, we neglect 
the existence of islands of explodability. 

What happens at intermediate metallicity between solar and zero? Predicting 
what happens to a metal-free star is relatively simple, because its evolution does 
not depend on the interplay between metals and stellar winds. The fate of a solar 
metallicity star is more problematic, because we must account for line-driven stellar 
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winds, but most observational data about stellar winds are for nearly solar metallic- 
ity stars and allow us to calibrate our models for such high metallicity. In contrast, 
modelling intermediate metallicities is significantly more complicated, because the 
details depend on the interplay between metals and stellar winds and only limited 
data are available for calibration (mostly data for the Large and Small Magellanic 
Clouds, e.g. [307]). 

As a rule of thumb (see e.g. [133] 322]), we can draw the following considera- 
tions. If the zero-age main sequence (ZAMS) mass of a star is large (mzams < 30 
Mo), then the amount of mass lost by stellar winds is the main effect which deter- 
mines the mass of the compact remnant. At low metallicity (< 0.1 Zo) and for a 
low Eddington factor (Iz < 0.6), mass loss by stellar winds is not particularly large. 
Thus, the final mass mg, and the carbon-oxygen mass mco of the star may be suf- 
ficiently large to avoid a core-collapse SN explosion: the star may form a massive 
BH ( Z 20 Mo) by direct collapse, unless a pair-instability or a pulsational-pair in- 
stability SN occurs. At high metallicity (~ Zo) or large Eddington factor (Iz > 0.6), 
mass loss by stellar winds is particularly efficient and may lead to a small mg, and 
mco: the star is expected to undergo a core-collapse SN and to leave a relatively 
small remnant. 

If the ZAMS mass of a star is relatively low (8 < mzams < 30 Mo), then stellar 
winds are not important (with the exception of asymptotic giant branch stars), re- 
gardless of the metallicity. In this case, the details of the SN explosion (e.g. energy 
of the explosion and amount of fallback) are crucial to determine the final mass of 
the remnant. This general sketch may be affected by several factors, such as pair- 
instability SNe, pulsational pair-instability SNe (e.g. [375]) and an island scenario 
for core-collapse SNe (e.g. [110]). 

The effect of pair-instability and pulsational pair-instability SNe is clearly shown 
in Figure[1.3] The top panel was obtained accounting only for stellar evolution and 
core-collapse SNe. In contrast, the bottom panel also includes pair-instability and 
pulsational pair-instability SNe. This figure shows that the mass of the compact 
remnant strongly depends on the metallicity of the progenitor star if mzams < 30 
Mo. In most cases, the lower the metallicity of the progenitor is, the larger the max- 
imum mass of the compact remnant [163}|225]|48] [232] 233] [322] [321]. However, for 
metal-poor stars (Z « 1073) with ZAMS mass 230 > mzams/Mo > 110 pair insta- 
bility SNe lead to the complete disruption of the star and no compact remnant is left. 
Only very massive (mzams > 230 Mo) metal-poor (Z < 10-2) stars can collapse to 
a BH directly, producing intermediate-mass BHs (i.e, BHs with mass Z 100 Mo). 

If Z < 10? and 110 > mzams 2 60 Mo, the star enters the pulsational pair- 
instability SN regime: mass loss is enhanced and the final BH mass is smaller 
(mpg ~ 30 — 55 Mo, bottom panel of Fig.|1.3) than we would have expected from 
direct collapse (mgg ~ 50 — 100 Mo, top panel of Fig. [1.3]. 

Finally, the mass spectrum of relatively low-mass stars (8 < mzams < 30 Mc) 
is not significantly affected by metallicity. The assumed core-collapse SN model is 
the most important factor in this mass range [133]. 

The models presented here do not take into account stellar rotation. Studied in 
[239], stellar rotation produces larger cores because of chemical mixing. 
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This shifts the minimum ZAMS mass for a star to undergo pulsational pair instabil- 
ity and PISN to lower values, because of the larger He core masses. The net result 
is a decrease of the maximum BH mass for fast-rotating stars (7 150 km s7!) with 
respect to low-rotating stars [236]. 
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Fig. 1.3 Mass of the compact remnant (M;em) as a function of the ZAMS mass of the star (mzams). 
Lower (upper) panel: pulsational pair-instability and pair-instability SNe are (are not) included. In 
both panels: dash-dotted brown line: Z = 2.0 x 10-2; dotted dark orange line: Z = 1.7 x 1072; 
dashed red line: Z = 1.4 x 102; solid red line: Z = 1.0 x 1072; short dash-dotted orange line: 
Z = 8.0x 102; short dotted light orange line: Z = 6.0 x 1073; short dashed green line: Z = 4.0 x 
10-3; dash-double dotted green line: Z = 2.0 x 1073; dash-dotted light blue line: Z = 1.0 x 1073; 
dotted blue line: Z = 5.0 x 1074; dashed violet line: Z = 2.0 x 1074. A delayed core-collapse SN 
mechanism has been assumed, following the prescriptions of [133]. This Figure was adapted from 
Figures 1 and 2 of Spera & Mapelli (2017, [321]). 
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1.2.5 Compact object spins 


In the previous Sections, we have seen that the connection between the mass of a 
BH and the properties of its progenitor star is still highly uncertain. Our knowledge 
on the origin of BH spins is even more uncertain. It is reasonable to assume that 
a compact object inherits the spin of its progenitor (or at least of the core of its 
progenitor) if the progenitor collapses to a BH directly, without any SN explosion. 
In contrast, mass ejection during a SN explosion can significantly dissipate part of 
the final spin of the progenitor star. Hence, if the final spin of the progenitor star is 
not negligible, we would expect large birth spins for the most massive BHs, which 
form from direct collapse, and low birth spins for NSs and light BHs, which form 
from successful SN explosions. 

Observations of Galactic pulsars seem to confirm the idea that NSs are born with 
relatively low x. Young pulsars (i.e., pulsars which did not have much time to slow 
down or to be recycled) have a value of the spin parameter y < 0.01 [185][219]. This 
indicates that most of the spin of the progenitor was lost either before or during the 
SN. In contrast, millisecond pulsars are significantly spun up by mass accretion from 
a companion star: the fastest millisecond pulsar (B1937+21, [166]) has y < 0.4, but 
this value has almost nothing to do with the initial spin of the NS. 

On the other hand, the observed spins of BHs make us "scratch our head". Rel- 
atively low-mass BHs in high-mass X-ray binaries seem to show very large spins 
[249]. Considering that these BHs did not have enough time 
to spin up by mass accretion (for the short lifetime of their companion stars), this 
suggests large spins at birth. In contrast, LIGO- Virgo results support low values of 
Xett for most of the BBHs [7], which can be interpreted either as strongly misaligned 
spins with respect to the binary angular momentum or as very low values of 7. 

From a theoretical perspective, the key question is then: what is the final spin of 
the progenitors of BHs and NSs? The final spin of a massive single star depends on 
mass loss and angular momentum transport. A mildly efficient transport by merid- 
ional currents (as adopted in, e.g., [106] [203]), leads to a non-negligible final spin 
of the star and to a high spin of the BH, if born from direct collapse [52] [236]. In 
contrast, [136] investigate efficient angular momentum transport via the magnetic 
Tayler-Spruit instability and find extremely slow spins (x ~ 10-2) for 
BHs born from single stars. 

If the magnetic Tayler-Spruit instability is the dominant process for angular mo- 
mentum transfer in massive stars, BH spins larger than y ~ 107? can be produced 
only by tidal torques in binary stars or by chemically homoge- 
neous evolution [382]. Overall, the initial rotation speed of very massive 
stars and the process of angular momentum transfer remain uncertain, hampering 
the predictive power of theoretical models on the spin distribution of BHs. 
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1.2.6 Natal kicks 


Compact objects are expected to receive a natal kick from the parent SN explosion, 
because of asymmetries in the neutrino flux and/or in the ejecta (see fora 
review). The natal kick has a crucial effect on the evolution of a binary compact 
object, because it can either unbind the binary or change its orbital properties. For 
example, a SN kick can increase the orbital eccentricity or misalign the spins of the 
two members of the binary. 

Unfortunately, it is extremely difficult to quantify natal kicks from state-of-the- 
art SN simulations and measurements of natal kicks are scanty, especially for BHs. 
As to NSs, indirect observational estimates of SN kicks give contrasting results. 
Hobbs et al. (2005, (169]) found that a single Maxwellian with root mean square 
Occsn = 265 km s^! can match the proper motions of 73 young single pulsars 
in the Milky Way. Other works suggest a bimodal velocity distribution, with a first 
peak at low velocities (e.g. ~ 0 km s^! according to or ~ 90 km s^! according 
to [34]) and a second peak at high velocities (> 600 km s~! according to or 
~ 500 km s^! for (34]). Similarly, Verbunt et al. (2017, (357]) indicate that a double 
Maxwellian distribution provides a significantly better fit to the observed velocity 
distribution than a single Maxwellian. Finally, the analysis of Beniamini & Piran 
(2016, [55]) shows that low kick velocities (< 30 km s7!) are required to match the 
majority of Galactic BNSs, especially those with low eccentricity. 

A possible interpretation of these observational results is that natal kicks depend 
on the SN mechanism (e.g. electron-capture versus core-collapse SN, e.g. [146]) or 
on the binarity of the NS progenitor. For example, if the NS progenitor evolves in 
a close binary system (i.e., in a binary system where the two stars have exchanged 
mass with each other, see Section[I.3.1) it might undergo an ultra-stripped SN (see 
[341] and references therein for more details). A star can undergo an ultra-stripped 
SN explosion only if it was heavily stripped by mass transfer to a companion 
[343]. The natal kick of an ultra-stripped SN should be low [343], because of the 
small mass of the ejecta (< 0.1 Mc). Low kicks (< 50 km s7!) for ultra-stripped 
core-collapse SNe are also confirmed by recent hydrodynamical simulations 
177]. 

As to BHs, the only indirect measurements of natal kicks arise from spatial distri- 
butions, proper motions and orbital properties of BHs in X-ray binaries (e.g. [250]). 
Evidence for a relatively small natal kick has been found for both GRO J1655-40 
[372] and Cygnus X-1 [374], whereas H 1705-250 and XTE J1118 +480 
require high kicks (> 100 km s^ !). By analysing the position of BHs in 
X-ray binaries with respect to the Galactic plane, Repetto et al. (2012, (291]) sug- 
gest that BH natal kicks should be as high as NS kicks. Repetto et al. (2017, [292]) 
perform a similar analysis but accounting also for binary evolution, and find that at 
least two BHs in X-ray binaries (H 1705-250 and XTE J1118 +480) require high 
kicks. 

Most models of BBH evolution assume that natal kicks of BHs are drawn from 
the same distribution as NS kicks, but reduced by some factor. For example, linear 
momentum conservation suggests that 
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m 
PNG (1.7) 
MBH 


where vgy is the natal kick of a BH with mass mpg, (mys) © 1.33Mo is the average 
mass of NSs and vys is randomly drawn from the NS kick distribution (e.g., 
from a single or double Maxwellian). 

Alternatively, the natal kick can be reduced by the amount of fallback, under the 
reasonable assumption that fallback quenches the initial asymmetries. Following 
[133], 

vBH = (1 — fø) vys, (1.8) 


where ff, quantifies the fallback (ff, = 0 for no fallback and fp = 1 for direct 
collapse). Most studies assume that BHs born from direct collapse receive no kick 
[133]. 

The model by can unify BH kicks and NS kicks, naturally accounting for 
ultra-stripped SNe and electron-capture SNe, which mostly lead to reduced kicks 
(e.g., [335] BA1]). According to this toy model, the kick of a compact object can be 
described as 
(mws) Mej 


Mrem (mej ) i 


vk = fos (1.9) 
where em is the mass of the considered compact remnant, Mej is the mass of the 
ejecta, (mj) is the ejecta mass of a SN that leaves a single NS with mass (mys), 
and fgos is a number randomly drawn from a Maxwellian probability density curve 
with one-dimensional root-mean square velocity dispersion o1p = 265 km s^! 
[169]. To derive this model, assume that the Maxwellian distribution fitted 
by is a good proxy to the distribution of the kicks of NSs born from single 
stars and that kicks are the effect of asymmetries in mass ejecta (e« Mej, see also 
[65] [66] [337]), modulated by linear momentum conservation (e mzzl). This simple 
formalism seems to solve the tension between the local merger rate of BNSs derived 
from gravitational-wave interferometers [7] [6] and the proper motions of Galactic 
young pulsars [169]. 

In summary, natal kicks are one of the most debated issues about compact ob- 
jects. Their actual amount has dramatic implications on the merger rate and on the 
properties (spin and mass distribution) of merging compact objects. 


1.3 Binaries of stellar black holes 


Naively, one could think that if two massive stars are members of a binary system, 
they will eventually become a BBH and the mass of each BH will be the same 
as if its progenitor star was a single star. This is true only if the binary system 
is sufficiently wide (detached binary) for its entire evolution. If the binary is tight 
enough, it will evolve through several processes which might significantly change 
its final fate. 
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The so-called binary population-synthesis codes have been used to investigate the 
effect of binary evolution processes on the formation of BBHs in isolated binaries 
(e.g. 
[340]. These are semi-analytic codes which combine a description 
of stellar evolution with prescriptions for SN explosions and with a formalism for 
binary evolution processes. In the following, we mention some of the most important 
binary-evolution processes and we briefly discuss their treatment in the most used 
population-synthesis codes. 


1.3.1 Mass transfer 


If two stars exchange matter to each other, it means they undergo a mass transfer 
episode. This might be driven either by stellar winds or by an episode of Roche- 
lobe filling. When a massive star loses mass by stellar winds, its companion might 
be able to capture some of this mass. This will depend on the amount of mass which 
is lost and on the relative velocity of the wind with respect to the companion star. 
Based on the Bondi & Hoyle (1944, [61]) formalism, Hurley et al. (2002, [172]) 
describe the mean mass accretion rate by stellar winds as 


1 Gm; NV? Ow 1 
mo = mil, 1.10 
: V1—e2 ( ve ) 2a? [1 T (Vorb/Vw)?]3/2 | | i ! 


where e is the binary eccentricity, G is the gravitational constant, m2 is the mass of 
the accreting star, vw is the velocity of the wind, «w ~ 3/2 is an efficiency constant, 
a is the semi-major axis of the binary, vog = \/G (mı 4- m2)/a is the orbital velocity 
of the binary (m, being the mass of the donor), and m; is the mass loss rate by the 
donor. Since |7; | is usually quite low (|r;| < 107? Mo yr^!) and vy is usually quite 
high (> 1000 km s^! for a line-driven wind) with respect to the orbital velocity, this 
kind of mass transfer is usually rather inefficient. However, most of the observed 
high-mass X-ray binaries and in particular all the systems with a WR 
star companion are wind-fed systems. 

Mass transfer by Roche lobe overflow is usually more efficient than wind accre- 
tion. The Roche lobe of a star in a binary system is a teardrop-shaped equipotential 
surface surrounding the star. The Roche lobes of the two members of the binary 
are connected in just one point, which is the Lagrangian L1 point. A widely used 
approximate formula for the Roche lobe is 


0.49 42/3 
a 
0.64255 + In (1-- 41/3)" 


Rui = (1.11) 


where a is the semi-major axis of the binary and q = mı /m (m, and m are the 
masses of the two stars in the binary). This formula describes the Roche lobe of a 
star with mass mı, while the corresponding Roche lobe of a star with mass mz (Ry?) 
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is obtained by swapping the subscripts. A star overfills (underfills) its Roche lobe 
when its radius is larger (smaller) than the Roche lobe. If a star overfills its Roche 
lobe, a part of its mass flows toward the companion star which can accrete (a part 
of) it. 

Mass transfer obviously changes the mass of the two stars in a binary, and thus 
the final mass of the compact remnants of such stars, but also the orbital properties 
of the binary. If mass transfer is non conservative, which is the most realistic case, 
it leads to an angular momentum loss, which in turn affects the semi-major axis. 
Recently, [63] show that the hypothesis of a highly non-conservative mass transfer 
(mass accretion efficiency fur < 0.5) is in tension with LVC data, if we assume that 
all BBHs observed by the LVC form via isolated binary evolution. 

A crucial information about Roche lobe overflow is whether it is stable or unsta- 
ble and on which timescale. The most commonly used approach can be described 
as follows [369] [283] [348] [172] [105]. Let us assume that the stellar radius and mass 
are connected by a simple relation R « m5 . Thus, the variation of the donor’s radius 
during Roche lobe is 

Ri dm, 
d — ot = d m; dt ` 


(1.12) 


In the above equation, the term 2m is due to nuclear burning, while the term with € 
measures the adiabatic or thermal response of the donor star to mass loss. Note that 
mi is the mass loss from the donor; hence it is always negative. 
Similarly, the change of the size of the Roche lobe of the donor RL, can be 
estimated as 
dR, 4 oR, ı Rui dm, 


dt ot ` 6L mj dt’ 


(1.13) 


where in depends on tides and GW radiation, while jj, describes the response of 
the Roche lobe to mass loss: the Roche lobe might shrink or expand. If & > C, the 
Roche lobe shrinks faster than the radius of the star does and the mass transfer is 
unstable, otherwise it remains stable until the radius changes significantly by nuclear 
burning. 

Mass transfer can be unstable either on a dynamical timescale (if € describes the 
adiabatic response of the donor and 6 < 61) or on a thermal timescale (if 6 describes 
the thermal response of the donor and ¢ < C1). If mass transfer is dynamically 
unstable or both stars overfill their Roche lobe, then the binary is expected to merge 
— if the donor lacks a steep density gradient between the core and the envelope -, or 
to enter common envelope (CE) - if the donor has a clear distinction between core 
and envelope. 


1.3.2 Common envelope (CE) 


If two stars enter CE, their envelope(s) stop co-rotating with their cores. The two 
stellar cores (or the compact object and the core of the companion star, if the binary 
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is already single degenerate) are embedded in the same non-corotating envelope and 
start spiralling in as an effect of gas drag exerted by the envelope. Part of the orbital 
energy lost by the cores as an effect of this drag is likely converted into heating of 
the envelope, making it more loosely bound. If this process leads to the ejection of 
the envelope, then the binary survives, but the post-CE binary is composed of two 
naked stellar cores (or of a compact object and a naked stellar core). Moreover, the 
orbital separation of the two cores (or the orbital separation of the compact object 
and the stellar core) is considerably smaller than the initial orbital separation before 
the CE, as an effect of the spiral ip] This circumstance is crucial for the fate of a 
BBH. In fact, if the binary which survives a CE phase evolves into a BBH, this BBH 
will have a short semi-major axis (a S 100 Re), much shorter than the sum of the 
maximum radii of the progenitor stars, and may be able to merge by GW emission 
within a Hubble time. 

In contrast, if the envelope is not ejected, the two cores (or the compact object 
and the core) spiral in till they eventually merge. This premature merger of a binary 
during a CE phase prevents the binary from evolving into a BBH. The cartoon in 
Figure[1.4]summarizes these possible outcomes. 

The œ formalism is the most common formalism adopted to describe a 
common envelope. The basic idea of this formalism is that the energy needed to 
unbind the envelope comes uniquely from the loss of orbital energy of the two cores 
during the spiral in. The fraction of the orbital energy of the two cores which goes 
into unbinding the envelope can be expressed as 


G 1 1 
AE = a (Ept — Epi) =e Me? ( ) , (1.14) 
2 af di 


where Ep ; (Ept) is the orbital binding energy of the two cores before (after) the CE 

phase, aj (ar) is the semi-major axis before (after) the CE phase, m; and mc» are the 

masses of the two cores, and œ is a dimensionless parameter that measures which 

fraction of the removed orbital energy is transferred to the envelope. If the primary 

is already a compact object (as in Figure[1.4). Mc2 is the mass of the compact object. 
The binding energy of the envelope is 


G Meny, M1 | Meny,2 M2 
À Ry R5 i 


Eenv = (1.15) 
where mı and m are the masses of the primary and the secondary member of the 
binary, Meny,1 and Meny,2 are the masses of the envelope of the primary and the 
secondary member of the binary, R, and Rz are the radii of the primary and the 
secondary member of the binary, and À is the parameter (or the function) which 


? A short-period (from a few hours to a few days) binary system composed of a naked helium core 
and BH might be observed as an X-ray binary, typically a WR X-ray binary. In the local Universe, 
we know a few (c 7) WR X-ray binaries, in which a compact object (BH or NS) accretes mass 
through the wind of the naked stellar companion (see e.g. for more details). These rare X-ray 
binaries are thought to be good progenitors of merging compact-object binaries. 
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Fig. 1.4 Schematic representation of the evolution of a BBH through CE. The companion of the 
BH is initially in the main sequence (MS). In the cartoon, the BH is indicated by a black circle, 
while the MS companion is indicated by the light blue circle. When the companion evolves off the 
MS, becoming a giant star, it overfills its Roche lobe. The BH and the giant star enter a CE (the CE 
is indicated in orange, while the core of the giant is represented by the dark blue circle). The core 
of the giant and the BH spiral in because of the gas drag exerted by the envelope. If the envelope is 
ejected, we are left with a new binary, composed of the BH and the naked helium core of the giant. 
The new binary has a much smaller orbital separation than the initial binary. If the naked helium 
core becomes a BH and its natal kick does not disrupt the binary, then a BBH is born, possibly 
with a small semi-major axis. In contrast, if the envelope is not ejected, the BH and the helium 
core spiral in, till they merge together. A single BH is left in this case. 


measures the concentration of the envelope (the smaller A is, the more concentrated 
is the envelope). 

By imposing AE = Ee, we can derive the value of the final semi-major axis 
ag for which the envelope is ejected. This means that the larger (smaller) œ is, the 
larger (smaller) the final orbital separation. If the resulting a; is lower than the sum 
of the radii of the two cores (or than the sum of the Roche lobe radii of the cores), 
then the binary will merge during CE, otherwise the binary survives. 

Actually, we have known for a long time (see for a review) that this simple 
formalism is a poor description of the physics of CE, which is considerably more 
complicated. A healthy treatment of CE should take into account not only the or- 
bital energy of the cores and the binding energy of the envelope, but also i) the 
thermal energy of the envelope, which is the sum of radiation energy and Kinetic 
energy of gas particles [161], ii) the recombination energy (as the envelope expands 
it cools down, the plasma recombines and some atoms even form molecules, releas- 
ing binding energy, (189]), iii) tidal heating/cooling from stellar spin down/up [175], 
iv) nuclear fusion energy [176], v) the enthalpy of the envelope [174], and vi) the 
accretion energy, which might drive outflows and jets [319]|213][212][320] 207] [89]. 


1 Formation channels of single and binary stellar-mass black holes 21 


Moreover, the envelope concentration parameter À cannot be the same for all 
stars. It is expected to vary wildly not only from star to star but also during different 
evolutionary stages of the same star. Several authors have estimated Eeny 
directly from their stellar models, significantly improving this formalism. However, 
even in this case, we cannot get rid of the œ parameter. 

Thus, it is essential to model the physics of CE with analytic models and nu- 
merical simulations. A lot of effort has been put on this in the last few years, but 
the problem remains largely unconquered. Several recent studies investigate the on- 
set of CE, when an unstable mass transfer prevents the envelope from co-rotating 
with the core and leads to the plunge-in of the companion inside the envelope 
[210] [211) [208} [209] [358]. 

Several hydrodynamical simulations model the fast spiral in phase after plunge- 
in [263], when the two cores spiral in on a dynamical time scale 
(~ 100 days). At the end of this dynamical spiral in only a small fraction of the 
envelope (~ 25%, (263]) appears to be ejected in most simulations. When the two 
cores are sufficiently close that they are separated only by a small gas mass, the 
spiral in slows down and the system evolves on the Kelvin-Helmholtz timescale of 
the envelope (~ 10?—? years). 

Simulating the system for a Kelvin-Helmholtz timescale is prohibitive for cur- 
rent three-dimensional simulations (e.g., [198]). Fragos et al. (2019, (129]) reduce 
the complexity of the problem by simulating the entire CE evolution in just one di- 
mension, with the hydrodynamic stellar evolution code MESA [2702711272]. They 
evolve a binary system composed of a NS and a 12 Mo red super-giant star for a 
thermal timescale. In their model, envelope ejection is mostly driven by the thermal 
energy of the envelope and by the orbital energy, while recombination contributes 
to < 10% of the total energy required to eject the envelope. The the final system is 
a NS - naked helium star system with an orbital separation of a few Ro, i.e. a good 
progenitor for a merging BNS. Their simulated system can be reproduced by the a 
formalism if œ ~ 5. The apparently unphysical value of œ > 1 is motivated by the 
fact that orbital energy is only one of the energy sources that participate in envelope 
ejection. 


1.3.3 Alternative evolution to CE 


Massive fast rotating stars can have a chemically homogeneous evolution (CHE): 
they do not develop a chemical composition gradient because of the mixing induced 
by rotation. This is particularly true if the star is metal poor, because stellar winds are 
not efficient in removing angular momentum. If a binary is very tight, the spins of 
its members are even increased during stellar life, because of tidal synchronisation. 
The radii of stars following CHE are usually much smaller than the radii of stars 
developing a chemical composition gradient [220]. This implies that even very 
tight binaries (few tens of solar radii) can avoid CE. 
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Marchant et al. (2016, (238]) simulate tight binaries whose components are fast 
rotating massive stars. A number of their simulated binaries evolve into contact 
binaries where both binary components fill and even overfill their Roche volumes. 
If metallicity is sufficiently low and rotation sufficiently fast, these binaries may 
evolve as "over-contact" binaries: the over-contact phase differs from a classical CE 
phase because co-rotation can, in principle, be maintained as long as material does 
not overflow the L2 point. This means that spiral-in can be avoided, resulting in a 
stable system evolving on a nuclear timescale. 

Such over-contact binaries maintain relatively small stellar radii during their evo- 
lution (few ten solar radii) and may evolve into a BBH with a very short orbital 
period. This scenario predicts the formation of merging BHs with relatively large 
masses (> 20 Mo), nearly equal mass (q = 1), and with large aligned spins (i.e., 
large Yes). While a large positive Yer is not consistent with most of the LIGO- 
Virgo BBHs, the CHE model is a viable alternative to explain systems with large 
and aligned spins. Moreover, mass loss during the naked helium core phase might 
reduce the progenitor's spin and lead to lower BH spins than predicted by early 


models [296]. 


1.3.4 BBH spins in the isolated binary evolution model 


Most evolutionary processes in an isolated binary star (tides, mass transfer, CE) lead 
to the alignment of the spins of the components to the orbital angular momentum of 
the binary [172]. The only? evolutionary process that can significantly misalign BH 
spins with respect to the orbital angular momentum is the SN explosion 
[299] [143] [2]. 

Under the assumption that the system has no time to re-align spins between the 
first and the second SN and that SNe do not change the direction of compact-object 
spins but only the direction of the orbital angular momentum of the binary, we can 
derive the angle between the direction of the spins of the two compact objects and 
that of the orbital angular momentum of the binary system as 


cos 0 = cos (V1) cos(v2) +sin(Vv1) sin (v2) cos (@), (1.16) 


where v; (with i = 1, 2) is the angle between the new (L5) and the old (L) orbital 
angular momentum after a SN (i — 1 corresponding to the first SN, i — 2 corre- 
sponding to the second SN), so that cos (v) = Lnew - Egg, while @ is the phase of 
the projection of Ê into the orbital plane. As shown by [299], the most commonly 
adopted SN kick models fail producing a significant misalignment. 

Hence, we expect that the isolated binary evolution model has a preference for 
BH spins aligned with the orbital angular momentum of the binary, especially in the 
case of massive BHs which undergo just a failed SN. 


3 Stegmann & Antonini (2021, (328]) recently proposed a possible spin flip mechanism during 
mass transfer. 
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1.3.5 Summary of the isolated binary formation channel 


In this Section, we have highlighted the most important aspects and the open issues 
of the isolated binary formation scenario, i.e. the model which predicts the forma- 
tion of merging BHs through the evolution of isolated binaries. For isolated binaries 
we mean stellar binary systems which are not perturbed by other stars or compact 
objects. 

To summarize, let us illustrate schematically the evolution of an isolated stel- 
lar binary (see e.g. [329] [147]) which can give birth to merging BHs like 
GW150914 and the other massive BHs observed by the LVC [51/4] [3] [10] [8] [7]. In 
the following discussion, several details of stellar evolution have been simplified or 
skipped to facilitate the reading for non specialists. 

Figure shows the evolution of an isolated binary system composed of two 
massive stars. These stars are gravitationally bound since their birth. Initially, the 
two stars are both on the main sequence (MS). When the most massive one leaves 
the MS, its radius starts inflating and can grow by a factor of a hundreds. The most 
massive star becomes a giant star with a helium core and a large hydrogen enve- 
lope. If its radius equals the Roche lobe (equation [11]. the system starts a stable 
mass-transfer episode. Some mass is lost by the system, some is transferred to the 
companion. After several additional evolutionary stages, the primary collapses to a 
BH. A direct collapse is preferred with respect to a SN explosion if we want the BH 
mass to be > 20 M; (as most BBHs observed by the LVC) and if we want the binary 
to remain bound (direct collapse implies almost no kick, apart from neutrino mass 
loss). At this stage, the semi-major axis of the system is still quite large: hundreds 
to thousands of solar radii. 

When the secondary also leaves the MS, growing in radius, the system enters a 
CE phase: the BH and the helium core spiral in. If the orbital energy and the thermal 
energy of the envelope are not sufficient to unbind the envelope, then the BH merges 
with the helium core leaving a single BH. In contrast, if the envelope is ejected, we 
are left with a new binary system, composed of the BH and of a stripped naked 
helium star. The new binary system has a much smaller orbital separation (tens of 
solar radii) than the pre-CE binary, because of the spiral in. If this new binary system 
remains bound after the naked helium star undergoes a SN explosion or if the naked 
helium star is sufficiently massive to directly collapse to a BH, the system evolves 
into a tight BBH, which might merge within a Hubble time. 

The most critical quantities in this scenario are the masses of the two stars and 
also their initial orbital separation (with respect to the stellar radii): a BBH can 
merge within a Hubble time only if its initial orbital separation is tremendously 
small (few tens of solar radii, unless the eccentricity is rather extreme, [274]), but 
a massive star (> 20 Mo) can reach a radius of several thousand solar radii during 
its evolution. Thus, if the initial orbital separation of the binary star is tens of solar 
radii, the binary system merges before it can become a BBH. On the other hand, if 
the initial orbital separation is too large, the two BHs will never merge. In this sce- 
nario, the two BHs can merge only if the initial orbital separation of the progenitor 
binary star is in the range which allows the binary to enter CE and then to leave 
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Fig. 1.5 Schematic evolution of an isolated binary star which can give birth to a BBH merger (see 


e.g. [50][228][329][147]). 
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a short period BBH. This range of initial orbital separations dramatically depends 
on CE efficiency and on the details of stellar mass and radius evolution. Two pos- 
sible alternatives to CE are very efficient stable mass transfer and CHE 
[238] 220]. 

There is one more interesting comment to add. Some stellar wind models predict 
a maximum BH mass of ~ 65 — 70 Mo from single stellar evolution and from the 
collapse of the residual hydrogen envelope at solar metallicity (e.g., [145]). 
This means that the maximum total mass of a BBH is mpop = mı +m ~ 130 — 140 
Mo. However, the maximum total mass of a BBH which merges within a Hubble 
time is only mror ~ 80 — 90 Mo (145], as shown in Figure[I.6] The reason for this 
difference is that if the initial orbital separation is sufficiently small (a ~ 10? — 104 
Ro, (323]) to allow mass transfer and CE, then the two stars lose almost all their 
hydrogen envelope: they might become a BBH that merges within a Hubble time, 
but their total mass is at most equal to the total mass of the two naked helium cores. 
In contrast, if the two stars are metal poor (Z < 0.0002) and if the initial orbital 
separation is too large to induce mass transfer and CE, the binary star becomes a 
BBH almost without mass loss, with a total mass up to 140 Mo, but the final semi- 
major axis is too large to allow coalescence by GW emission. 


1.4 The dynamics of binary black holes (BBHs) 


In the previous Sections of this review, we discussed the formation of BBHs as 
isolated binaries. There is an alternative channel for BBH formation: the dynamical 
evolution scenario. 


1.4.1 Dynamically active environments 


Collisional dynamics is important for the evolution of binaries only if they are in 
a dense environment ( Z 10? stars pe?) such as a star cluster. On the other hand, 
astrophysicists believe that the vast majority of massive stars (which are BH pro- 
genitors) form in star clusters ([196][370] [371] [282], but see also (367]). 

There are several different flavours of star clusters. Globular clusters are 
old stellar systems (~ 12 Gyr), mostly very massive (total mass Msc > 10* Ma) 
and dense (central density pe > 10^ Mo pc~’). They are sites of intense dynamical 
processes, such as the gravothermal catastrophe. Globular clusters represent a small 
fraction of the baryonic mass in the local Universe ( 5 1 per cent, ([162]). Most 
studies of dynamical formation of BBHs focus on globular clusters (e.g. 
[224] [08] 99] [54] [338] [67] [305] [300] [299] 6] 304] [03] [35] [30] [124] [298] B01 [186]). 

Young star clusters are young ( < 100 Myr), relatively dense (px 10? Mo pc?) 
stellar systems, and are thought to be the most common birthplace of massive stars 
[282]. When they disperse by gas evaporation or are disrupted by the tidal 
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Fig. 1.6 Upper panel: total mass (mror = m + m2) of all BBHs formed from a stellar population 
of 2 x 107 massive binary stars simulated with MOBSE and with different metallicities (from 0.0002 
to 0.02) . Lower panel: total mass of the BBHs merging within a Hubble time from the same 


simulations. See https: //mobse-webpage.netlify.app/|for more details on MOBSE 


and on these runs. Courtesy of Nicola Giacobbo. 


field of their host galaxy, their stellar content is released into the field. Thus, it 
is reasonable to expect that a large fraction of BBHs which are now in the field 
may have formed in young star clusters, where they participated in the dynamics 
of the cluster. A fraction of young star clusters might survive gas evaporation and 
tidal disruption and evolve into older open clusters, like M67. Figure shows 
a snapshot of an N-body simulation of a young star cluster. Population-synthesis 
prescriptions are included in this simulation to follow the evolution of binary stars 
and the formation of BBHs. Studies of BBHs in young and open clusters include 
(2801 40] [386] 223] [181] [58] [39] [135] [193] 194] |4 1] |93] [94] [95]. 

Finally, nuclear star clusters lie in the nuclei of galaxies. Nuclear star clusters 
are rather common in galaxies (e.g. [258]), are usually more mas- 
sive and denser than globular clusters, and may co-exist with super-massive BHs 
(SMBHs). Stellar-mass BHs formed in the innermost regions of a galaxy could even 
be "trapped" in the accretion disc of the central SMBH, triggering their merger (see 
e.g. [244]. These features make nuclear star clusters unique among star 
clusters, for the effects that we will describe in the next Sections. 
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Fig. 1.7 Snapshot of an N-body simulation of a young star cluster. Different particle sizes and 
colors refer to different star luminosities and temperatures, as estimated by population-synthesis 
calculations. This simulation has been included in [233]. 


1.4.2 Two-body encounters, dynamical friction and core-collapse 


The main driver of the dynamics of star clusters is gravity force. Gravitational two- 
body encounters between stars lead to local fluctuations in the potential of the star 
cluster and drive major changes in the internal structure of the star cluster over a 
two-body relaxation timescale [325] [326]: 


3 


oO 
trix = 0.34 —— 
ic Beda pA 


(1.17) 


where o is the local velocity dispersion of the star cluster, (m) is the average stel- 
lar mass in the star cluster, p is the local mass density, G is the gravity force and 
InA ~ 10 is the Coulomb logarithm. The two-body relaxation timescale is the time 
needed for a typical star in the stellar system to completely lose memory of its initial 
velocity due to two-body encounters. In star clusters, fj, is much shorter than the 
Hubble time (tj ~ 10— 100 Myr in young star clusters, (282]), while in galaxies 
and large-scale structures it is much longer than the lifetime of the Universe. Hence, 
close encounters are common in dense star clusters. 
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Dynamical friction is another consequence of gravity force: a massive body of 
mass M orbiting in a sea of lighter particles feels a drag force that slows down its 
motion over a timescale [74]: 


3 o? 


4 (22)? Gina Mp(r) 


tpr(M) = (1.18) 


It is apparent that two-body relaxation and dynamical friction are driven by the 
same force and are related by 


tpe(M) ~ s (1.19) 


i.e., dynamical friction happens over a much shorter timescale than two-body re- 
laxation and leads to mass segregation (or mass stratification) in a star cluster. This 
process speeds up the collapse of the core of a star cluster and can trigger the so- 
called Spitzer's instability [324]. Two-body relaxation, dynamical friction and their 
effects play a crucial role in shaping the demography of BBHs in star clusters, as 
we discuss below. 


1.4.3 Binary — single encounters 


We now review what are the main dynamical effects which can affect a BBH, start- 
ing from binary — single star encounters. Binaries have a energy reservoir, their 


internal energy: 

mee a (1.20) 

2 r 

where u = m, m»/ (m 4- m3) is the reduced mass of the binary (whose components 
have mass mı and m2), v is the relative velocity between the two members of the 
binary, and r is the distance between the two members of the binary. As shown by 
Kepler's laws, Eim = —Ep = —G mı m / (2a), where Ey is the binding energy of the 
binary system, a being its semi-major axis. 

The internal energy of a binary can be exchanged with other stars only if the 
binary undergoes a close encounter with a star, so that its orbital parameters are 
perturbed by the intruder. This happens only if a single star approaches the binary 
by few times its orbital separation. We define this close encounter between a binary 
and a single star as a three-body encounter. For this to happen with a non-negligible 
frequency, the binary must be in a dense environment, because the rate of three- 
body encounters scales with the local density of stars. Three-body encounters have 
crucial effects on BH binaries, such as hardening, exchanges and ejections. 
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1.4.4 Hardening 


If a BBH undergoes a number of three-body encounters during its life, we expect 
that its semi-major axis will shrink as an effect of the encounters. This process is 
called dynamical hardening. 

Following [165], we call hard binaries (soft binaries) those binaries with binding 
energy larger (smaller) than the average kinetic energy of a star in the star cluster. 
According to Heggie's law [165], hard binaries tend to harden (i.e., to become more 
and more bound) via binary-single encounters. In other words, a fraction of the 
internal energy of a hard binary can be transferred into kinetic energy of the intruder 
and of the centre-of-mass of the binary during three body encounters. This means 
that the binary loses internal energy and its semi-major axis shrinks. 

Most BBHs are expected to be hard binaries, because BHs are among the most 
massive bodies in star clusters. Thus, BBHs are expected to harden as a consequence 
of three-body encounters. The hardening process may be sufficiently effective to 
shrink a BBH till it enters the regime where GW emission is efficient: a BBH which 
is initially too loose to merge may then become a GW source thanks to dynami- 
cal hardening. The hardening rate for hard binaries with semi-major axis a can be 


estimated as [165] 
d /1 p 
—|-)=27G¢— 1.21 
a (+) 22645. (121) 
where é ~ 0.1 — 10 is a dimensionless hardening parameter (which has been esti- 
mated through numerical experiments, [167] ), p is the local mass density of 
stars, © is the local velocity dispersion and G is the gravity constant. 

Dynamical hardening is the main responsible for the shrinking of a binary, till 


it reaches a semi-major axis short enough for GW emission to become effective, 
which can be derived with the following equation [274]: 


da 64 G^mym»; (m +m) 3 
= : 1.22 
d 5 gcü-eyme ^ i 


By combining equations and it is possible to make a simple analytic 
estimate of the evolution of the semi-major axis of a BBH which is affected by 
three-body encounters and by GW emission: 


da 
dt 


64 G? m, m (mj +m) "E 
5. gGü-eyn , 


ipae n (1.23) 
oO 
This equation holds under the assumption that the binary star is hard, the total mass 
of the binary star is much larger than the average mass of a star in the star cluster 
(exchanges are unlikely) and that most three-body encounters have a small impact 
parameter. The first part of the right-hand term of equation [1.23] accounts for the 
effect of three-body hardening on the semi-major axis. It scales as da/dt « — a?, 
indicating that the larger the binary is, the more effective the hardening. This can 
be easily understood considering that the geometric cross section for three body 
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interactions with a binary scales as a?. The second part of the right-hand term of 
equation accounts for energy loss by GW emission. It is the first-order ap- 
proximation of the calculation by Peters (1964, [274]). It scales as da /dt« — a 
indicating that GW emission becomes efficient only when the two BHs are very 
close to each other. 
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Fig. 1.8 Time evolution of the semi-major axis of three BH binaries estimated from equation|1.23 
Blue dashed line: BH binary with masses mı = 200 Mo, m = 30 Mo; red solid line: m; = 36 
Mo, m = 29 Mo; black dot-dashed line: m; = 14 Mo, m» = 7.5 Mo. For all BH binaries: § = 1, 
p = 10? Ma pc 2, o = 10 km s^l, e = 0 (here we assume that p, o and e do not change during 
the evolution), initial semi-major axis of the BH binary a; — 10 AU. 


In Figure[1.8] we solve equation [1.23] numerically for three BBHs with different 
mass. All binaries evolve through i) a first phase in which hardening by three body 
encounters dominates the evolution of the binary, ii) a second phase in which the 
semi-major axis stalls because three-body encounters become less efficient as the 
semi-major axis shrink, but the binary is still too large for GW emission to become 
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efficient, and iii) a third phase in which the semi-major axis drops because the binary 
enters the regime where GW emission is efficient. 


1.4.5 Exchanges 


Dynamical exchanges are three-body encounters during which one of the members 
of the binary is replaced by the intruder. Exchanges may lead to the formation of new 
BBHs: if a binary composed of a BH and a low-mass star undergoes an exchange 
with a single BH, this leads to the formation of a new BBH. This is a fundamental 
difference between BHs in the field and in star clusters: a BH which forms as a 
single object in the field has negligible chances to become member of a binary 
system, while a single BH in the core of a star cluster has good chances of becoming 
member of a binary by exchanges. 

Exchanges are expected to lead to the formation of many more BBHs than they 
can destroy, because the probability for an intruder to replace one of the members 
of a binary is ~ 0 if the intruder is less massive than both binary members, while it 
suddenly jumps to ~ 1 if the intruder is more massive than one of the members of 
the binary [168]. Since BHs are among the most massive bodies in a star cluster after 
their formation, they are very efficient in acquiring companions through dynamical 
exchanges. Thus, exchanges are a crucial mechanism to form BH binaries dynam- 
ically. By means of direct N-body simulations, Ziosi et al. (2014, ([386]) show that 
> 90 per cent of BBHs in young star clusters form by dynamical exchange. More- 
over, BBHs formed via dynamical exchange will have some distinctive features with 
respect to field BBHs (see e.g. [386]): 


* BBHs formed by exchanges will be (on average) more massive than isolated 
BBHs, because more massive intruders have higher chances to acquire compan- 
ions; 

* exchanges trigger the formation of highly eccentric BBHs; eccentricity is then 
significantly reduced by circularisation induced by GW emission, if the binary 
enters the regime where GW emission is effective; 

* BBHs born by exchange will likely have misaligned spins: exchanges and other 
dynamical interactions tend to lead to isotropically distributed spin directions 
with respect to the binary orbital plane, because dynamical interactions remove 
any memory of previous alignments. 


Zevin et al. (2017, ([384]) compare a set of simulations of field binaries with a set 
of simulations of globular cluster binaries, run with the same population-synthesis 
code. The most striking difference between merging BH binaries in their globular 
cluster simulations and in their population-synthesis simulations is the dearth of 
merging BHs with mass < 10 Mo in the globular cluster simulations. This is likely 
due to the fact that exchanges tend to destroy binaries composed of light BHs. 

Figure[1.9] compares the masses of simulated BBH mergers in isolated binaries 
and in young star clusters. While the maximum total mass of BBH mergers in iso- 
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Fig. 1.9 Mass of the secondary BH versus mass of the primary BH in a set of BBH mergers. 
Yellow stars: O1 and O2 BBHs plus GW190412. Blue, green and red symbols: simulated BBHs 
in young star clusters with metallicity Z — 0.0002, 0.002 and 0.02. Filled symbols: dynamical 
exchanges. Open symbols: original binaries. Gray contours: BBHs from isolated binary evolution. 
Simulations from . Courtesy of Ugo N. Di Carlo. 


lated binaries is ~ 80 Mo (see Figure[I.6), BBH mergers in young star clusters have 
total masses up to ~ 130 Mo. This difference springs mainly from two facts: 


1) single stellar evolution at low metallicity can lead to the formation of single BHs 
with mass up to ~ 70 Mo, if the hydrogen envelope collapses to the final BH. 
In the field such massive BHs remain single, while in a star cluster they can 
acquire companions via dynamical exchanges and they can merge by dynamical 
hardening and GW emission; 
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ii) in star clusters, stellar collisions are quite common and can lead to the build 
up of more massive BHs, even with mass inside the pair instability mass gap 


290]. 


Spin misalignments are another possible feature to discriminate between field 
binaries and star cluster binaries (e.g. [115]. We expect that an isolated bi- 
nary system in which both the primary and the secondary component undergo direct 
collapse results in a BBH with nearly aligned spins. For dynamically formed BH 
binaries we expect misaligned, or even nearly isotropic spins, because any original 
spin alignment is completely reset by three-body encounters. 

Currently, we have limited constraints on BH spins from GW detections. In a 
few events, such as GW151226 [4], GW170729 [8] and especially GW 190412 [17], 
the measured value of x. is significantly larger than zero, indicating at least partial 
alignment. The vast majority of events listed in GWTC-2 have y. consistent with 
zero. This might be the result of either low spins or misaligned spins or a combi- 
nation of both. On the other hand, the most recent population study by the LVC 
[15] shows that ~ 12% to 44% of BBH systems have spins tilted by more than 90° 
with respect to their orbital angular momentum, supporting a negative effective spin 
parameter. 

Interestingly, GW 190412 has a non-zero precessing spin %p to the 90% credible 
level and GW190521 [20], the most massive BBH event to date, shows mild 
evidence for non-zero precessing spin. This might support a dynamical formation 


for these events (e.g., ). 


1.4.6 Stellar mergers and BHs in the pair-instability mass gap 


Binary-single encounters in young star clusters trigger collisions between massive 
stars. This process is sped up by dynamical friction, which can generate an excess 
of massive stars in the cluster core in less than a Myr. Stellar collisions can lead 
to the formation of very massive stars (> 150 Mo), blue straggler stars and 
unusually massive BHs. 

Figure[I. 10]shows the dynamical formation of a BBH with primary mass m = 88 
Mo and secondary mass m» = 47.5 Mo [95]. The masses of the components of this 
BBH are similar to those of GW190521 (18][20]. In particular, the primary BH has 
a mass in the pair-instability mass gap. Its formation is possible because its stellar 
progenitor is the result of the merger between a giant star with a well-developed 
He core and a main sequence companion (MS). The result of this stellar merger is 
a massive core helium burning (cHeB) star with an over-sized hydrogen envelope 
with respect to the He core [93] [94] [95] [187] [290]. Given the short timescale for He, 
C, O, Ne and Si burning with respect to H burning, the star collapses to a BH before 
the He core grows above the threshold for pair instability: the result is a ~ 88 Mo 
BH. After its formation, this BH acquires a companion by dynamical exchanges and 
merges within a Hubble time. 
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Fig. 1.10 Cartoon of the dynamical assembly of a GW190521-like BBH from the simulations of 
. Courtesy of Ugo N. Di Carlo. 


This is a viable scenario to produce systems like GW190521 not only because 
of the good match of the BBH mass, but also because a BH born from the direct 
collapse of a very massive star might form with a large spin (depending on the 
final spin of the massive progenitor) and because this dynamical formation leads 
to isotropically oriented spins. The main uncertainties of this model are mass loss 
during stellar collision and the possible loss of a fraction of the envelope 
during the final collapse [118]. 
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1.4.7 Direct three-body binary formation 


In the most massive star clusters (globular clusters, nuclear star clusters), stellar 
velocities are so large that dynamical encounters can unbind most of the original 
binary stars (i.e., those binary stars that were already there at the formation of the 
star cluster). The minimum relative velocity v; between a binary star and an intruder 
star needed to unbind a binary is in fact [317]: 


ass "Em mo (mi +m im) (1:24) 


ma (mj +m2)a 


where m, m2 and m3 are the masses of the two binary members and that of the 
intruder, while a is the binary semi-major axis. 

In these extreme environments, most BBHs are expected to form by direct en- 
counters of three single bodies [304], during core collapse. This leads to the 
formation of extremely hard BBHs, which survive further ionization from intruders. 
The timescale for binary formation via three single body encounters is [199] [25]: 


fa = 0.1 Myr [ 2 q Z i mau Y^ (1.25) 
pcc VYE A T06 pe 30kms!) (30M) ^ ' 


where n is the local stellar density, o is the local velocity dispersion of the star clus- 
ter and mpy is the typical BH mass in the star cluster. The properties of BBHs born 
from three single body encounters are similar to those of BBHs born via dynamical 
exchanges: they tend to be more massive than isolated binaries, have high initial 
eccentricity and isotropically oriented spins (e.g. (27]). 

Direct three-body encounters are likely the most common BBH formation chan- 
nel in globular clusters and nuclear star clusters [253], while binary — single star 
exchanges are likely the most common formation channel of BBHs in young star 
clusters [93]. In young star clusters, dynamical exchanges affect both already 
formed BHs and their stellar progenitors. 


1.4.8 Dynamical ejections 


During three-body encounters, a fraction of the internal energy of a hard binary 
is transferred into kinetic energy of the intruders and of the centre-of-mass of the 
binary. As a consequence, the binary and the intruder recoil. The recoil velocity is 
generally of the order of a few km s~!, but can be up to several hundred km s~!. 
Since the escape velocity from a globular cluster is ~ 30 km s^! and the escape 
velocity from a young star cluster or an open cluster is even lower, both the recoiling 
binary and the intruder can be ejected from the parent star cluster. If the binary 
and/or the intruder are ejected, they become field objects and cannot participate in 
the dynamics of the star cluster anymore. Thus, not only the ejected BBH stops 
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hardening, but also the ejected intruder, if it is another compact object, loses any 
chance of entering a new binary by dynamical exchange. 

A general expression for the recoil velocity of the binary center of mass if (m; + 
mo) >> (m) (where (m) is the average mass of a star in a star cluster) is 


EG) Vi 28 Fh, (1.26) 


mı +m m, +m + (m)) 


where Ey = Gm, m2/(2a) is the binary binding energy. The above equation can 
help us deriving the minimum binding energy above which a binary star is ejected 
by a binary-single encounter £y, min [246]: 


(mı +m)’ 2 


DE (m)? Vesc: (1.27) 


Ep, min n 


where vesc is the escape velocity from the star cluster. A BBH will merge inside the 
star cluster only if E min > Ey, Gw, where Ey gw is the minimum binding energy to 
reach coalescence by GW emission. 

Most BNSs, BBHs and BH-NS systems in young star clusters are ejected be- 
fore they merge [93] [289]. Dynamical ejections of BNS and BH-NS binaries were 
proposed to be one of the possible explanations for the host-less short gamma-ray 
bursts, i.e. gamma-ray bursts whose position in the sky appears to be outside any 
observed galaxy [122]. Host-less bursts may be ~ 25% all short gamma-ray bursts. 

In general, ejections of compact objects and compact-object binaries from their 
parent star cluster can be the result of at least three different processes: 


* dynamical ejections (as described above); 
e SN kicks [169/1133]; 
* GW recoil [205] 151]. 


GW recoil is a relativistic kick occurring when a BBH merges. It is the result 
of asymmetric linear momentum loss by GW emission, when the binary has asym- 
metric component masses and/or misaligned spins. It results in kick velocities up to 
thousands of km s~! and usually of the order of hundreds of km s^ !. 

Ejections by dynamics, SN kicks or GW recoil may be the main process at work 
against mergers of second-generation BHs, where for second-generation BHs we 
mean BHs which were born from the merger of two BHs rather than from the col- 
lapse of a star [143]. In globular clusters, open clusters and young star clusters, a BH 
has good chances of being ejected by three-body encounters before it merges 
and a very high chance of being ejected by GW recoil after it merges [301][234]. The 
only place where merging BHs can easily avoid ejection by GW recoil are nuclear 
star clusters, whose escape velocity can reach hundreds of km s^! [25][28][127]234]. 
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1.4.9 Formation of intermediate-mass black holes by runaway 
collisions 


In Section [1.2.4] we have mentioned that intermediate-mass black holes (IMBHs, 
i.e. BHs with mass 100 < mgu < 104 Mo) might form from the direct collapse of 
metal-poor extremely massive stars [321]. Other formation channels have been pro- 
posed for IMBHs and most of them involve the dynamics of star clusters. The 
formation of massive BHs by runaway collisions has been originally proposed 
about half a century ago and was then elaborated by several authors (e.g. 
(279] 281) 278) [160] [130] 226] 231) [1491 223] [93] [95] 297]). 


Fig. 1.11 Cartoon of the runaway collision scenario in dense young star clusters (see e.g. [281]). 
From left to right: (1) the massive stars (red big stars) and the low-mass stars (yellow small stars) 
follow the same initial spatial distribution; (2) dynamical friction leads the massive stars to sink 
to the core of the cluster, where they start colliding between each other; (3) a very massive star 
C» 100 Mo) forms as a consequence of the runaway collisions; (4) this massive star might be able 
to directly collapse into a BH. 


The basic idea is the following (as summarized by the cartoon in Figure [I.T1). 
In a dense star cluster, dynamical friction makes massive stars to decelerate 
because of the drag exerted by lighter bodies, on a timescale tpr(M) ~ ((m) /M) tax 
(equation{I.18). 

Since the two-body relaxation timescale in a young star cluster can be as short 
as tix ~ 10 — 100 Mo [282], for a star with mass M > 40 Mo we estimate tpp < 
2.5 Myr: dynamical friction is very effective in dense massive young star clusters. 
Because of dynamical friction, massive stars segregate to the core of the cluster 
before they become BHs. 

If the most massive stars in a dense young star cluster sink to the centre of the 
cluster by dynamical friction on a time shorter than their lifetime (i.e., before core- 
collapse SNe take place, removing a large fraction of their mass), then the density 
of massive stars in the cluster core becomes extremely high. This makes collisions 
between massive stars extremely likely. Actually, direct N-body simulations show 
that collisions between massive stars proceed in a runaway sense, leading to the 
formation of a very massive (> 100 Mo) star [281]. The main open question is: 
“What is the final mass of the collision product? Is the collision product going to 
collapse to an IMBH?". 
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There are essentially two critical issues: i) how much mass is lost during the 
collisions? ii) how much mass does the very-massive star lose by stellar winds? 

Hydrodynamical simulations of colliding stars show that massive star 
can lose ~ 25% of their mass during collisions. Even if we optimistically assume 
that no mass is lost during and immediately after the collision (when the colli- 
sion product relaxes to a new equilibrium), the resulting very massive star will be 
strongly radiation-pressure dominated and is expected to lose a significant fraction 
of its mass by stellar winds. Recent studies including the effect of the Eddington 
factor on mass loss show that IMBHs cannot form from runaway colli- 
sions at solar metallicity. At lower metallicity (Z < 0.1 Zo) approximately 10 — 3096 
of runaway collision products in young dense star clusters can become IMBHs by 
direct collapse (they also avoid being disrupted by pair-instability SNe). 

The majority of runaway collision products do not become IMBHs but they end 
up as relatively massive BHs (~ 20 — 90 Mo, [223]). If they remain inside their par- 
ent star cluster, such massive BHs are extremely efficient in acquiring companions 
by dynamical exchanges. Mapelli (2016, ([223]) find that all stable binaries formed 
by the runaway collision product are BBHs and thus are possibly important sources 
of GWs in the LIGO- Virgo range. 


1.4.10 Hierarchical BBH formation and IMBHs 


The runaway collision scenario occurs only in the early stages of the evolution of a 
star cluster, when massive stars are still alive (the lifetime of a ~ 30 Mo star is ~ 6 
Myr). However, it has been proposed that IMBHs form even in old clusters (e.g., 
globular clusters) by repeated mergers of smaller BHs (e.g. 
[235]. 

The simple idea is illustrated in Figure [I-12] A stellar BBH in a star cluster is 
usually a rather hard binary. Thus, it shrinks by dynamical hardening till it may 
enter the regime where GW emission is effective. In this case, the BBH merges 
leaving a single more massive BH. Given its relatively large mass, the new BH has 
good chances to acquire a new companion by exchange. Then, the new BBH starts 
hardening again by three-body encounters and the story may repeat several times, 
till the main BH becomes an IMBH. 

This scenario has one big advantage: it does not depend on stellar evolution, 
so we are confident that the BH will grow in mass by mergers, if it remains in- 
side the cluster. However, there are several issues. First, the BBH may be ejected 
by dynamical recoil, received as an effect of three-body encounters. Recoils get 
stronger and stronger, as the orbital separation decreases (equation|1.27). The BBH 
will avoid ejection by dynamical recoil only if it is sufficiently massive ( Z 50 Me 
for a dense globular cluster, [85]). If the BBH is ejected, the loop breaks and no 
IMBH is formed. 

Second and even more important, the merger of two BHs involves a relativistic 
kick. This kick may be as large as hundreds of km s~! (205], leading to the ejection 
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Hardening induced by 3-body encounters 


Merger induced by GWs 


The merger product acquires a new companion through exchange 


Fig. 1.12 Cartoon of the repeated merger scenario in old star clusters (see, e.g., [246][149]). From 
top to bottom and from left to right: (1) a BBH undergoes three-body encounters in a star cluster; 
(2) three-body encounters harden the BBH, shrinking its semi-major axis; (3) the BBH hardens 
by three-body encounters till it enters the regime where GW emission is efficient: the BBH semi- 
major axis decays by GW emission and the binary merges; (4) a single bigger BH forms as result 
of the merger and may acquire a new companion by dynamical exchange (if it is not ejected by GW 
recoil); (5) the new BBH containing the bigger BH starts shrinking again by three-body encounters 
(1). This loop may be repeated several times till the main BH becomes an IMBH. 


of the BH from the parent star cluster [170]. Also in this case, the loop breaks and 
no IMBH is formed. Finally, even if the BH binary is not ejected, this scenario is 
relatively inefficient: if the seed BH is ~ 50 Mo, several Gyr are required to form 
an IMBH with mass ~ 500 Mo [246]. 

Monte Carlo simulations by Giersz et al. (2015, ([149]) show that both the run- 
away collision scenario and the repeated merger scenario can be at work in star 
clusters: runaway collision IMBHs form in the first few Myr of the life of a star 
cluster and grow in mass very efficiently, while repeated-merger IMBHs start form- 
ing much later ( Z 5 Gyr) and their growth is less efficient. 


1.4.11 Alternative models for massive BHs and IMBH formation in 
galactic nuclei 


Several additional models predict the formation of IMBHs in galactic nuclei. For 
example, Miller & Davies (2012, [247]) propose that IMBHs can efficiently grow 
in galactic nuclei from runaway tidal capture of stars, provided that the velocity 
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dispersion in the nuclear star cluster is 2 40 km s~!. Below this critical velocity, 
binary stars can support the system against core collapse, quenching the growth 
of the central density and leading to the ejection of the most massive BHs. Above 
this velocity threshold, the stellar density can grow sufficiently fast to enhance tidal 
captures and BH-star collisions. Tidal captures are more efficient than BH-star col- 
lisions in building up IMBHs, because the mass growth rate of the former scales 
as MIMBH & mimg“ (where mpypuy is the IMBH seed mass), whereas the mass 
growth of the latter scales as rit]upg ec mimpy [332]. 

Furthermore, McKernan et al. (2012, [243], see also [380] 381]) 
suggest that IMBHs could efficiently grow in the accretion disc of a SMBH. Nuclear 
star cluster members trapped in the accretion disk are subject to two competing ef- 
fects: orbital excitation due to dynamical heating by encounters with other stars and 
orbital damping due to gas drag. Gas damping is expected to be more effective than 
orbital excitation, quenching the relative velocity between nuclear cluster members 
and enhancing the collision rate. This favours the growth of IMBHs via both gas 
accretion and multiple stellar collisions. This mechanism might be considered as a 
gas-aided hierarchical merger scenario. 


As a final remark, it is worth mentioning that all the IMBH formation scenar- 
ios we have discussed here — i)runaway collisions of massive stars, ii) hierarchical 
merger of BHs and iii) BH trapping in AGN discs- can also work as possible forma- 
tion mechanisms for GW190521-like events: they can lead to the formation of BHs 
with mass in the pair instability mass gap and with large spins. These BHs are born 
in dense stellar environments, where they can acquire companions dynamically and 
form BBHs with isotropically oriented spins [20]. 


1.4.12 Kozai-Lidov resonance 


Unlike the other dynamical processes discussed so far, Kozai-Lidov (KL) resonance 
can occur both in the field and in star clusters. KL resonance appears 
whenever we have a stable hierarchical triple system (i.e., a triple composed of an 
inner binary and an outer body orbiting the inner binary), in which the orbital plane 
of the outer body is inclined with respect to the orbital plane of the inner binary. 
Periodic perturbations induced by the outer body on the inner binary cause 1) the 
eccentricity of the inner binary and ii) the inclination between the orbital plane of 
the inner binary and that of the outer body to oscillate. The semi-major axis of the 
binary star is not affected, because KL resonance does not imply an energy exchange 
between inner and outer binary. KL oscillations may enhance BBH mergers, because 
the timescale for merger by GW emission strongly depends on the eccentricity e of 
the binary taw œ (1 —e?)7/? (see equation[1.22] [274]). 

It might seem that hierarchical triples are rather exotic systems. This is not the 
case. In fact, ~ 10% of low-mass stars are in triple systems [345] [346] 287]. This 
fraction gradually increases for more massive stars [100], up to ~ 50% for B-type 
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stars [347]. In star clusters, stable hierarchical triple systems may form 
dynamically, via four-body or multiple-body encounters. 

Kimpson et al. (2016, (181]) find that KL resonance may enhance the BBH 
merger rate by ~ 40% in young star clusters and open clusters. According to Fra- 
gione et al. (2019, [125]), the merger fraction of BBHs in galactic nuclei can be up 
to ~ 5 — 8 times higher for triples than for binaries. On the other hand, Antonini et 
al. (2017, [26]) find that KL resonance in field triples can account for < 3 mergers 
Gpc 2? yr !. The main signature of the merger of a KL system is the non-zero eccen- 
tricity until very few seconds before the merger. Eccentricity might be significantly 
non-zero even when the system enters the LIGO- Virgo frequency range. 

KL resonances have an intriguing application in nuclear star clusters. If the stellar 
BH binary is gravitationally bound to the super-massive BH (SMBH) at the centre of 
the galaxy, then we have a peculiar triple system where the inner binary is composed 
of the stellar BH binary and the outer body is the SMBH [29]. Also in this case, the 
merging BBH has good chances of retaining a non-zero eccentricity till it emits 
GWs in the LIGO-Virgo frequency range. 


1.4.13 Summary of dynamics and open issues 


In this Section, we have seen that dynamics is a crucial ingredient to understand 
BBH demography. Dynamical interactions (three- and few-body close encounters) 
can favour the coalescence of BBHs through dynamical hardening. New BBHs can 
form via dynamical exchanges and via direct three-body encounters or GW cap- 
tures. Hierarchical mergers of BHs in dense star clusters and in AGN disks can trig- 
ger the formation of unusually massive BHs (like GW190521) and even IMBHs. 
Stellar mergers might also lead to the formation of unusually massive BHs (like 
GW190521) and IMBHs in dense star clusters. All these dynamical processes sug- 
gest a boost of the BBH merger rate in dynamically active environments. 

Overall, dynamically formed BBHs are expected to be more massive than BBHs 
from isolated binary evolution, with higher initial eccentricity and with misaligned 
spins. Also, KL resonances favour the coalescence of BBHs with extremely high 
eccentricity, even close to the last stable orbit. On the other hand, three-body en- 
counters might trigger the ejection of binary compact objects from their natal envi- 
ronment, inducing a significant displacement between the birth place of the binary 
and the location of its merger. 

Figure [1.13] summarizes one of the possible evolutionary pathways of merging 
BBHs which originate from dynamics (the variety of this formation channel is too 
large to account for all dynamical channels mentioned above in a single cartoon). As 
in the isolated binary case, we start from a binary star. In the dynamical scenario, it is 
not important that this binary evolves through Roche lobe or CE (although this may 
happen). After the primary has turned into a BH, the binary undergoes a dynamical 
exchange: the secondary is replaced by a massive BH and a new BBH forms. The 
new binary system is not ejected from the star cluster and undergoes further three- 
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g. 1.13 Schematic evolution of a merging BBH formed by dynamical exchange (see, e.g., 
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body encounters. As an effect of these three-body encounters, the binary system 
hardens enough to enter the regime in which GW emission is efficient: the BBH 
merges by GW decay. 

We expect dynamics to be important for BBHs also because massive stars (which 
are the progenitors of BHs) form preferentially in young star clusters [282], which 
are dynamically active places. It is reasonable to expect that most BHs participate 
in the dynamics of their parent star cluster before being ejected by dynamical recoil 
or relativistic kicks. 


1.5 BBHs in the cosmological context 


At design sensitivity, LIGO and Virgo will observe BBHs up to redshift z ~ 1. The 
farthest event observed to date, GW 190521, is at z ~ 0.8, when the Universe was 
only ^ 6.8 Gyr old. Next generation ground-based GW detectors, such as the Ein- 
stein Telescope [284], will observe BBH mergers up to redshift z Z, 10, when the 
Universe was 5 500 Myr old. Almost the entire Universe will be transparent to 
BBH mergers and it will be possible to investigate the cosmic evolution of stars 
and galaxies through the observations of BBH mergers. In other words, we will do 
cosmology through BBH mergers. 

Accounting for the cosmological context of BBH mergers might appear as a des- 
perate challenge, because of the humongous dynamical range: the orbital separa- 
tions of BBHs are of the order of tens of solar radii, while cosmic structures are 
several hundreds of Mpc. Several theoretical studies have addressed this challenge, 
adopting two different methodologies. 


1.5.1 Data-driven semi-analytic models 


Some authors (e.g. 
[310]) combine the outputs of population synthesis codes with analytic pre- 
scriptions. The main ingredients are the cosmic star formation rate density and the 
evolution of metallicity with redshift [214] 215] [81] [8S0] [79]. In some previous work 
(e.g. [96] [197]) a Press-Schechter like formalism is adopted, to include the mass of 
the host galaxy in the general picture. Lamberts et al. (2016, [[197]) even include a 
redshift-dependent description for the mass-metallicity relation (hereafter MZR), to 
account for the fact that the mass of a galaxy and its observed metallicity are deeply 
connected. The main advantage of this procedure is that the star formation rate and 
the metallicity evolution can be derived more straightforwardly from the data. The 
main drawback is that it is more difficult to trace the evolution of the host galaxy of 
the BBH, through its galaxy merger tree. 

In the semi-analytic description, the merger rate density evolution can be de- 


scribed as [309]: 
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Zz 4 Zmax 
R(2) = Í v) we) | " Q(Z.g(z,sZ)dZ|d/, (128) 


Zmin 

where dr(z)/dz = (14-2) H(z)-, with H(z) = Ho [(12-2)) Qu 4- Q4]. 
Furthermore, Zmin and Zmax are the minimum and maximum metallicity, y(z’) is 
the cosmic star formation rate (SFR) at redshift z’, F (7, z, Z) is the rate of compact 
binaries that form at redshift z’ from stars with metallicity Z and merge at redshift 
z, and 7(Z) is the merger efficiency, namely the ratio between the total number 
"A ror(Z) of compact binaries (formed from a coeval population) that merge within 
an Hubble time (ty, < 14 Gyr) and the total initial mass M,(Z) of the simulation 
with metallicity Z. The value of F (7, z, Z) and that of ](Z) can be calculated ei- 
ther from catalogues of population-synthesis simulations or from phenomenological 
models. 

This formalism yields an evolution of the merger rate density with redshift as 
shown in Fig. [1.14] From the upper panel, we see that the cosmic merger rate den- 
sity evolution is sensitive to the choice of the CE parameter œ. From the lower 
panel, it is apparent that the BBH merger rate density is tremendously affected by 
the metallicity evolution and by the observational uncertainties on metallicity evo- 
lution. Most models agree on this result (e.g., [310]). 
The reason for this trend is twofold. On the one hand, according to current mod- 
els of BBH formation, the merger efficiency n is 2—4 orders of magnitude higher 
at low metallicity (Z « 0.0002) than at high metallicity. Hence, the merger rate is 
extremely sensitive to the underlying metallicity evolution. On the other hand, the 
uncertainties on metallicity evolution from observational data are large, as discussed 
by, e.g., [214][218][81] [80] [309][79]. 

The trend of the cosmic merger rate in Fig. |1.14|is similar to the trend of the 
cosmic star formation rate density curve, modulated by both metallicity evolution 
and time delay] The peak of the BBH merger rate density curve (z ~ 3 — 4, de- 
pending on o) is at a higher redshift than the peak of the cosmic star formation 
rate density (z = 2). This happens because the lower the metallicity is, the higher 
is the efficiency of BBH mergers, and metal-poor stars are more common at high 
redshift than at low redshift. With this approach, it will also be possible to put con- 
straints, in a statistical sense, on the fraction of BBH mergers from the dynamical 
channel with respect to BBH mergers from the isolated evolutionary channel (e.g., 


[298] [62] [309] 373] [385]). 


1.5.2 Cosmological simulations 


The alternative approach to reconstruct the BBH merger history feeds the outputs of 
population-synthesis simulations into cosmological simulations [265| 
[33] (32) 31], through a Monte Carlo approach. This has the 


^ We define as time delay the time elapsed between the formation of the progenitor binary and the 
merger of the two BHs. 
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Fig. 1.14 Cosmic merger rate density of BBHs in the comoving frame (Zpgpgg) as a function of the 
look-back time (fp, lower x—axis) and of the redshift (z, upper x—axis) from some of the models 
presented in [310], which are based on equation [1.28] We use catalogues of population synthesis 
simulations run with MOBSE [228] [147]. In both panels: thin gray line and right y— axis: cosmic 
star formation rate density as a function of redshift from [215]. Gray shaded box: 90% credible 
interval for the local merger rate density of BBHs, as inferred from the LVC data, considering the 
union of the rates obtained with model A, B and C in [7]. The width of the grey shaded area on the 
x—axis corresponds to the instrumental horizon obtained by assuming BBHs with total mass of 20 
Mo and O2 sensitivity [9]. Upper panel: BBH merger rate density for different assumptions on the 
a parameter (i.e. the efficiency of CE ejection). Models for œ = 0.5, 1, 2, 3, 5, 7 and 10 are shown. 
Lower panel: the solid red line shows the median value of the BBH merger rate density for ot = 5, 
while the shaded areas show the 50% credible areas estimated from the uncertainties on metallicity 
evolution (hatched region) and star formation rate density evolution (shaded region). The cosmic 
star formation rate density evolution is modelled as in and is obtained by fitting data. The 
metallicity evolution uses the fitting formula by [90], adapted as decribed in [309]. Courtesy of 
Filippo Santoliquido. 
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clear advantage that the properties of the host galaxies can be easily reconstructed 
across cosmic time. However, the ideal thing would be to have a high-resolution 
cosmological simulation (sufficient to resolve also small dwarf galaxies) with a box 
as large as the instrumental horizon of the GW detectors. This is obviously im- 
possible. High-resolution simulations have usually a box of a few comoving Mpc?, 
while simulations with a larger box cannot resolve dwarf galaxies. Moreover, this 
procedure requires to use the cosmic star formation rate density and the redshift- 
dependent MZR which are intrinsic to the cosmological simulations. While most 
state-of-the-art cosmological simulations reproduce the cosmic star formation rate 
density reasonably well, the MZR is an elusive feature, creating more than a trouble 
even in the most advanced cosmological simulations. 

For example, investigate the main properties of the host galaxies of merg- 
ing BHs, by combining their population-synthesis simulations with the 
ILLUSTRIS cosmological box (364] [365] [257]. The size of the ILLUSTRIS (length 
— 106.5 comoving Mpc) is sufficient to satisfy the cosmological principle, but 
galaxies with stellar mass < 108 Mo are heavily under-resolved. Their results show 
that BHs merging in the local Universe (z « 0.024) have formed in galaxies with rel- 
atively small stellar mass (< 10!° Mọ) and relatively low metallicity (Z < 0.1 Zo). 
These BHs reach coalescence either in the galaxy where they formed or in larger 
galaxies (with stellar mass up to ~ 10! Mo). In fact, most BHs reaching coales- 
cence in the local Universe appear to have formed in the early Universe ( Z 9 Gyr 
ago), when metal-poor galaxies were more common. A significant fraction of these 
high-redshift metal-poor galaxies merged within larger galaxies before the BBHs 
reached coalescence by GWs. Moreover, these models show that the mass spectrum 
and the other main properties of BBHs do not evolve significantly with redshift 
[229]. 

Schneider et al. (2017, (314], see also [159]) adopt a complementary ap- 
proach to study the importance of dwarf galaxies for GW detections. They use 
the GAMESH pipeline to produce a high-resolution simulation of the Local Group 
(length = 4 Mpc comoving). This means that the considered portion of the Universe 
is strongly biased, but the resolution is sufficient to investigate BH binaries in small 
( Z 10° Mo) dwarf galaxies. One of their main conclusions is that GW 150914-like 
events originate mostly from small metal-poor galaxies. 

Similarly, Cao et al. (2018, (70]) investigate the host galaxies of BBHs by ap- 
plying a semi-analytic model to the Millennium-II N-body simulation [64]. The 
Millennium-II N-body simulation is a large-box (length = 137 comoving Mpc) dark- 
matter only simulation. The physics of baryons is implemented through a semi- 
analytic model. Using a dark-matter only simulation coupled with a semi-analytic 
approach offers the possibility of improving the resolution significantly, but baryons 
are added only in post-processing. With this approach, Cao et al. (2018, (70]) find 
that BBHs merging at redshift z 5 0.3 are located mostly in massive galaxies (stellar 
mass Z2 x 10! Mo). 

Finally, combine population-synthesis simulations with the EAGLE 
cosmological simulation [312]. They show that the merger rate per galaxy strongly 
correlates with the stellar mass of a galaxy for both BBHs, BHNSs and BNSs. They 
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also find a dependence of the merger rate per galaxy on the star formation rate and 
a weaker dependence on the metallicity, but the correlation with the stellar mass is 
definitely the strongest one. This result has been used to optimize electromagnetic 
follow-up strategies, weighting the galaxies in the LVC uncertainty box by their stel- 
lar mass [101] [22]. These studies show that the combination of population-synthesis 
tools with cosmological simulations is an effective approach to understand the cos- 
mic evolution of the merger rate and the properties of the host galaxies of BBH 
mergers. 


1.6 Summary and outlook 


We reviewed our current understanding of the formation channels of BBHs. The last 
five years have witnessed considerable progress in this field, thanks to the ground- 
breaking results of the LVC and to a renewed interest in BH astrophysics, mostly 
triggered by GW data. 

The final stages of massive star evolution deserve particular attention to under- 
stand the mass and spins of BHs. GW data have recently challenged the concept of 
pair-instability mass gap, i.e. the existence of a mass range (between ~ 65 Mo and 
~ 120 Mo) in which we do not expect to find BHs as an effect of (pulsational) pair 
instability. 

The process of mass transfer between two massive stars in a binary system is 
another key aspect shaping the demography of BBHs. The efficiency of mass ac- 
cretion via Roche lobe overflow and the stability of mass transfer are possibly the 
main unknowns affecting the mass and the delay time of a BBH. In addition, the 
dynamical formation channel with its richness and multifaceted processes (harden- 
ing, exchanges, ejections, runaway collisions, hierarchical mergers, etc) adds to the 
complexity of the general picture. 

Current astrophysical models of BBH formation face a number of essential ques- 
tions: 


* impact of core overshooting, rotation and possibly magnetic fields on massive 
star evolution and BH formation; 

* angular momentum transfer in massive stars and its link to BH spins; 

e explodability of massive stars; 

* stability and efficiency of mass transfer; 

* evolution of common-envelope systems; 

* physics of stellar collisions and mergers; 

e impact of the environment (e.g., metallicity, star clusters versus field) on the for- 
mation of BBHs across cosmic time. 


Overall, we have a better understanding of the main dynamical processes with 
respect to both massive star evolution and SN physics. But our understanding of 
BBH dynamics is mainly hampered by two problems: i) we cannot properly model 
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the dynamics of BBHs if we do not know BH masses and spins from stellar evolu- 
tion; ii) dynamical simulations are still too computationally expensive to allow us to 
investigate the relevant parameter space. 

Next-generation ground-based detectors will observe BBH mergers up to redshift 
z © 10, beyond the epoch of cosmic reionization [179]. This will open a completely 
new scenario for the study of BBHs across cosmic time and for the characterization 
of their evolutionary pathways. 
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